module cam_history
!-------------------------------------------------------------------------------------------
! 
! The cam_history module provides the user interface for CAM's history output capabilities.
! It maintains the lists of fields that are written to each history file, and the associated
! metadata for those fields such as descriptive names, physical units, time axis properties,
! etc.  It also contains the programmer interface which provides routines that are called from
! the physics and dynamics initialization routines to define the fields that are produced by
! the model and are available for output, and the routine that is called from the corresponding
! run method to add the field values into a history buffer so that they may be output to disk.
!
! There are two special history files.  The initial file and the satellite track file.
!
! Public functions/subroutines:
!   addfld, add_default
!   intht
!   write_restart_history
!   read_restart_history
!   outfld
!   wshist
!   initialize_iop_history
!-----------------------------------------------------------------------

   use shr_kind_mod, only: r8 => shr_kind_r8, r4 => shr_kind_r4
   use spmd_utils,   only: masterproc, npes
   use ppgrid,       only: pcols
   use filenames, only: interpret_filename_spec
   use filenames, only: ncdata, bnd_topo
   use abortutils, only: endrun
   use pmgrid, only : dyndecomp_set

   use units,     only: getunit

   use hycoef, only: hyai, hybi, hyam, hybm, ps0

   use dyn_grid,     only: get_horiz_grid_dim_d, get_horiz_grid_d, get_dyn_grid_parm

   use pio,          only: file_desc_t, var_desc_t, io_desc_t, &
                           iotype_pnetcdf, iotype_netcdf, &
                           pio_noerr, pio_bcast_error, pio_internal_error, &
                           pio_seterrorhandling, pio_setdebuglevel, pio_setframe, &
                           pio_rearr_box, pio_rearr_none,  &
                           pio_nofill, pio_clobber,pio_offset, &
                           pio_int, pio_real, pio_double, pio_char, pio_offset, pio_unlimited, pio_global, &
                           pio_inq_dimid, pio_inq_dimlen, pio_inq_varid, pio_inq_varname, pio_inq_varndims, &
                           pio_def_dim, pio_def_var, pio_enddef, &
                           pio_put_att, pio_put_var, pio_get_att, pio_get_var, &
                           pio_write, pio_write_darray, &
                           pio_read_darray, pio_closefile, pio_freedecomp, pio_syncfile

   use scamMod,       only: scmlon,single_column
   use perf_mod,      only: t_startf, t_stopf
   use cam_logfile,   only: iulog	
   use cam_pio_utils, only: phys_decomp, dyn_decomp, dyn_stagger_decomp
   use cam_history_buffers, only : dim_index_3d,hbuf_accum_inst, hbuf_accum_add, hbuf_accum_min, hbuf_accum_max, &
        hbuf_accum_add00z, hbuf_accum_addlcltime
   use cam_history_support, only : max_fieldname_len, fieldname_suffix_len, max_chars, ptapes, fieldname_len, max_string_len, &
        date2yyyymmdd, registeredmdims, fieldname_lenp2, sec2hms, pflds, column_info, field_info, active_entry, hentry, &
        hist_mdims, register_hist_mdim, lookup_hist_mdim_indicies
   use sat_hist, only : is_satfile

   implicit none
   private

   public :: phys_decomp, dyn_decomp, dyn_stagger_decomp, fieldname_len ! forward from cam_pio_utils


!
! master_entry: elements of an entry in the master field list
!
   type master_entry
      type (field_info)     :: field            ! field information
      character(len=1)           :: avgflag(ptapes)  ! averaging flag
      character(len=max_chars) :: time_op(ptapes)  ! time operator (e.g. max, min, avg)
      logical               :: act_sometape     ! Field is active on some tape
      logical               :: actflag(ptapes)  ! Per tape active/inactive flag
      integer               :: htapeindx(ptapes)! This field's index on particular history tape
      type(master_entry), pointer :: next_entry ! The next master entry
   end type master_entry

   type (master_entry), pointer :: masterlinkedlist    ! master field linkedlist top

   type master_list
      type(master_entry), pointer :: thisentry
   end type master_list

   type (master_list), pointer :: masterlist(:) ! master field array for hash lookup of field

!  column_info structure has moved to cam_pio_utils



!


   type (active_entry), pointer :: tape(:)          ! history tapes
   type (active_entry), target,allocatable :: history_tape(:)          ! history tapes
   type (active_entry), target, allocatable :: restarthistory_tape(:)          ! restart history tapes


   type rvar_id
      type(var_desc_t), pointer :: vdesc
      integer :: type
      integer :: ndims
      integer :: dims(4)
      character(len=fieldname_lenp2) :: name
   end type rvar_id
   type rdim_id
      integer :: len
      integer :: dimid
      character(len=fieldname_lenp2) :: name      
   end type rdim_id
!
!   The size of these parameters should match the assignments in restart_vars_setnames and restart_dims_setnames below
!   
   integer, parameter :: restartvarcnt=29
   integer, parameter :: restartdimcnt=8
   type(rvar_id) :: restartvars(restartvarcnt)
   type(rdim_id) :: restartdims(restartdimcnt)



   integer :: nfmaster = 0             ! number of fields in master field list
   integer :: nflds(ptapes)            ! number of fields per tape

! per tape sampling frequency (0=monthly avg)

   integer :: i                        ! index for nhtfrq initialization
   integer :: nhtfrq(ptapes) = (/0, (-24, i=2,ptapes)/)  ! history write frequency (0 = monthly)
   integer :: mfilt(ptapes) = 30       ! number of time samples per tape
   integer :: nfils(ptapes)            ! Array of no. of files on current h-file
   integer :: ngroup(ptapes)           ! Array of no. of contiguous columns on current h-file
   integer :: nexcl(ptapes)            ! Actual number of excluded fields
   integer :: nincl(ptapes)            ! Actual number of included primary file fields
   integer :: ndens(ptapes) = 2        ! packing density (double (1) or real (2))
   integer :: ncprec(ptapes) = -999    ! netcdf packing parameter based on ndens
   real(r8) :: beg_time(ptapes)        ! time at beginning of an averaging interval
!
! Netcdf ids
!

   type(var_desc_t) :: gwid             ! var id for gaussian weights


   integer :: nscurf(ptapes)           ! First "current" second of day for each h-file
   integer :: ncsecf(ptapes)           ! First "current" second of date for each h-file

   logical :: rgnht(ptapes) = .false.  ! flag array indicating regeneration volumes
   logical :: hstwr(ptapes) = .false.  ! Flag for history writes
   logical :: empty_htapes  = .false.  ! Namelist flag indicates no default history fields
   logical :: htapes_defined = .false. ! flag indicates history contents have been defined

   integer :: luhrest = -1             ! history restart file unit number
   character(len=max_string_len) :: hrestpath(ptapes) = (/(' ',i=1,ptapes)/) ! Full history restart pathnames
   character(len=max_string_len) :: nfpath(ptapes) = (/(' ',i=1,ptapes)/) ! Array of first pathnames, for header
   character(len=max_string_len) :: cpath(ptapes)                   ! Array of current pathnames
   character(len=max_string_len) :: nhfil(ptapes)                   ! Array of current file names
   character(len=1)  :: avgflag_pertape(ptapes) = (/(' ',i=1,ptapes)/) ! per tape averaging flag
   character(len=8)  :: logname             ! user name
   character(len=16) :: host                ! host name
   character(len=max_string_len) :: ctitle = ' '      ! Case title
   character(len=8)  :: inithist = 'YEARLY' ! If set to '6-HOURLY, 'DAILY', 'MONTHLY' or
                                            ! 'YEARLY' then write IC file 
   character(len=fieldname_lenp2) :: fincl(pflds,ptapes) ! List of fields to add to primary h-file
   character(len=max_chars)       :: fincllonlat(pflds,ptapes) ! List of fields to add to primary h-file
   character(len=fieldname_lenp2) :: fexcl(pflds,ptapes) ! List of fields to rm from primary h-file
   character(len=fieldname_lenp2) :: fwrtpr(pflds,ptapes) ! List of fields to change default history output prec
   character(len=fieldname_suffix_len ) :: fieldname_suffix = '&IC' ! Suffix appended to field names for IC file

   character(len=22) ,parameter :: LT_DESC = 'mean (over local time)' ! local time description
   logical :: collect_column_output(ptapes)

   integer :: maxvarmdims=1
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  Hashing.
!
!  Accelerate outfld processing by using a hash function of the field name
!  to index masterlist and determine whehter the particular field is to
!  be written to any history tape.
!
!
!  Note: the outfld hashing logic will fail if any of the following are true:
!
!         1) The lower bound on the dimension of 'masterlist' is less than 1.
!
!         2) 'outfld' is called with field names that are not defined on
!            masterlist.  This applies to both initial/branch and restart
!            runs.
!
!         3) An inconsistency between a field's tape active flag
!            'masterlist(ff)%actflag(t)' and active fields read from
!            restart files.
!
!         4) Invoking function 'gen_hash_key' before the primary and secondary
!            hash tables have been created (routine bld_outfld_hash_tbls).
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!  User definable constants for hash and overflow tables.
!  Define size of primary hash table (specified as 2**size).
!
   integer, parameter :: tbl_hash_pri_sz_lg2 = 16
!
!  Define size of overflow hash table % of primary hash table.
!
   integer, parameter :: tbl_hash_oflow_percent = 20
!
!  Do *not* modify the parameters below.
!
   integer, parameter :: tbl_hash_pri_sz = 2**tbl_hash_pri_sz_lg2
   integer, parameter :: tbl_hash_oflow_sz = tbl_hash_pri_sz * (tbl_hash_oflow_percent/100.0_r8) 
!
!  The primary and overflow tables are organized to mimimize space (read:
!  try to maximimze cache line usage).
!
!  gen_hash_key(fieldname) will return an index on the interval
!  [0 ... tbl_hash_pri_sz-1].
!
!
!  Primary:
!  gen_hash_key(fieldname)-------+     +----------+
!                                |     |   -ii    | 1 ------>tbl_hash_oflow(ii)
!                                |     +----------+
!                                +-->  |    ff    | 2 ------>masterlist(ff)
!                                      +----------+
!                                      |          | ...
!                                      +----------+
!                                      |          | tbl_hash_pri_sz
!                                      +----------+
!
!  Overflow (if tbl_hash_pri() < 0):
!  tbl_hash_pri(gen_hash_key(fieldname))
!                         |
!                         |            +----------+
!                         |            |     1    | 1  (one entry on O.F. chain)
!                         |            +----------+
!                         |            |    ff_m  | 2
!                         |            +----------+
!                         +--------->  |     3    | 3  (three entries on chain)
!                                      +----------+
!                                      |    ff_x  | 4
!                                      +----------+
!                                      |    ff_y  | 5
!                                      +----------+
!                                      |    ff_z  | 6
!                                      +----------+
!                                      |          | ...
!                                      +----------+
!                                      |          | tbl_hash_oflow_sz
!                                      +----------+
!
!
   integer, dimension(0:tbl_hash_pri_sz-1) :: tbl_hash_pri ! Primary hash table
   integer, dimension(tbl_hash_oflow_sz) :: tbl_hash_oflow ! Overflow hash table
!
!  Constants used in hashing function gen_hash_key.
!  Note: if the constants in table 'tbl_gen_hash_key' below are modified,
!        changes are required to routine 'gen_hash_key' because of specific
!        logic in the routine that optimizes character strings of length 8.
!

   integer, parameter :: gen_hash_key_offset = z'000053db'

   integer, parameter :: tbl_max_idx = 15  ! 2**N - 1
   integer, dimension(0:tbl_max_idx) :: tbl_gen_hash_key = &
   (/61,59,53,47,43,41,37,31,29,23,17,13,11,7,3,1/)

   ! Used in output interpolation of unstructured grids
   logical :: interpolate, ncoloutput


!
!
! Generic procedures
!
!
! Public entities
!

!
! Filename specifiers for history, initial files and restart history files
! (%c = caseid, $y = year, $m = month, $d = day, $s = seconds in day, %t = tape number)
!
   character(len=max_string_len) :: rhfilename_spec = '%c.cam2.rh%t.%y-%m-%d-%s.nc' ! history restart
   character(len=max_string_len), public :: hfilename_spec(ptapes) = (/ (' ', i=1, ptapes) /) ! filename specifyer

! needed by history_default
   public :: init_masterlinkedlist

! Needed by initext

   public :: nhtfrq, mfilt, inithist, ctitle

! Needed by read_namelist

   public :: fincl, fincllonlat, fexcl, fwrtpr
   public :: pflds, ptapes, empty_htapes, ndens
   public :: avgflag_pertape, collect_column_output
   integer, public :: lcltod_start(ptapes) ! start time of day for local time averaging (sec)
   integer, public :: lcltod_stop(ptapes)  ! stop time of day for local time averaging, stop > start is wrap around (sec)

! Needed by stepon

   public :: hstwr
   public :: nfils

! Functions
   public :: init_restart_history     ! Write restart history data
   public :: write_restart_history     ! Write restart history data
   public :: read_restart_history      ! Read restart history data
   public :: wshist                    ! Write files out
   public :: outfld                    ! Output a field
   public :: intht                     ! Initialization
   public :: wrapup                    ! process history files at end of run
   public :: write_inithist            ! logical flag to allow dump of IC history buffer to IC file
   public :: addfld                    ! Add a field to history file
   public :: add_default               ! Add the default fields
   public :: get_hfilepath             ! Return history filename
   public :: get_ptapes                ! Return the number of tapes being used
   public :: get_hist_restart_filepath ! Return the full filepath to the history restart file
   public :: hist_fld_active           ! Determine if a field is active on any history file

CONTAINS

  subroutine init_masterlinkedlist()
    use dyn_grid, only : get_dyn_grid_parm
    integer :: plev, plevp

    nullify(masterlinkedlist)
    nullify(tape)
    plev = get_dyn_grid_parm('plev')
    plevp = get_dyn_grid_parm('plevp')

    call register_hist_mdim('lev',plev)
    call register_hist_mdim('ilev',plevp)

  end subroutine init_masterlinkedlist


   subroutine intht ()
!
!----------------------------------------------------------------------- 
! 
! Purpose: Initialize history file handler for initial or continuation run.
!          For example, on an initial run, this routine initializes "ptapes"
!          history files.  On a restart or regeneration  run, this routine 
!          only initializes history files declared beyond what existed on the 
!          previous run.  Files which already existed on the previous run have 
!          already been initialized (i.e. named and opened) in routine RESTRT.
! 
! Method: Loop over tapes and fields per tape setting appropriate variables and
!         calling appropriate routines
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
      use ioFileMod
      use shr_sys_mod, only: shr_sys_getenv
      use time_manager, only: get_prev_time, get_curr_time
      use cam_control_mod, only : nsrest
      use dycore,           only: dycore_is
      use sat_hist,        only: sat_hist_init
#if (defined SPMD)
      use spmd_utils, only : mpichar, mpicom
#endif
!
!-----------------------------------------------------------------------
!
! Local workspace
!
      integer :: t, f              ! tape, field indices
      integer :: begdim1           ! on-node dim1 start index
      integer :: enddim1           ! on-node dim1 end index
      integer :: begdim2           ! on-node dim2 start index
      integer :: enddim2           ! on-node dim2 end index
      integer :: begdim3           ! on-node chunk or lat start index
      integer :: enddim3           ! on-node chunk or lat end index
      type (dim_index_3d) :: dimind  ! 3-D dimension index
      integer :: day, sec          ! day and seconds from base date
      integer :: i, ii, ff         ! index
      integer :: rcode             ! shr_sys_getenv return code
      type(master_entry), pointer :: listentry
      character(len=32) :: fldname ! temp variable used to produce a left justified field name
                                   ! in the formatted logfile output

!
! Print master field list
!

      if (masterproc) then
         write(iulog,*) ' '
         write(iulog,*)' ******* MASTER FIELD LIST *******'
      end if
      listentry=>masterlinkedlist
      f=0
      do while(associated(listentry))
         f=f+1
         if(masterproc) then
            fldname = listentry%field%name
            write(iulog,9000) f, fldname, listentry%field%units, listentry%field%numlev, &
                             listentry%avgflag(1), listentry%field%long_name
9000        format(i5, 1x, a32, 1x, a16, 1x, i4, 1x, a1, 2x, 256a)
         end if
         listentry=>listentry%next_entry
      end do
      nfmaster = f
      if(masterproc) write(iulog,*)'BLDFLD:nfmaster=',nfmaster

!
!  Now that masterlinkedlist is defined and we are performing a restart run
!  (htapes_defined == .true.), construct primary and secondary hashing tables.
!
      if ( htapes_defined ) then
         call bld_outfld_hash_tbls()
         call bld_htapefld_indices()
         return
      end if
!
! Get users logname and machine hostname
!
      if ( masterproc )then
         logname = ' '
         call shr_sys_getenv ('LOGNAME',logname,rcode)
         host = ' '
         call shr_sys_getenv ('HOST',host,rcode)
      end if
#ifdef SPMD
! PIO requires netcdf attributes have consistant values on all tasks
      call mpibcast(logname, len(logname), mpichar, 0, mpicom)
      call mpibcast(host, len(host), mpichar, 0, mpicom)
#endif
!
! Override averaging flag for all fields on a particular tape if namelist input so specifies
!
      do t=1,ptapes
         if (avgflag_pertape(t) /= ' ') then
            call h_override (t)
         end if
      end do
!
! Define field list information for all history files.  
!
      call fldlst ()
!
! Loop over max. no. of history files permitted  
!
      if ( nsrest .eq. 3 ) then
         call get_prev_time(day, sec)  ! elapased time since reference date
      else
         call get_curr_time(day, sec)  ! elapased time since reference date
      end if
      do t=1,ptapes
         nfils(t) = 0            ! no. of time samples in hist. file no. t

! Time at beginning of current averaging interval.
         
         beg_time(t) = day + sec/86400._r8
      end do

!
! Initialize history variables
!
      do t=1,ptapes
         do f=1,nflds(t)
            begdim1  = tape(t)%hlist(f)%field%begdim1
            enddim1  = tape(t)%hlist(f)%field%enddim1
            begdim2  = tape(t)%hlist(f)%field%begdim2
            enddim2  = tape(t)%hlist(f)%field%enddim2
            begdim3  = tape(t)%hlist(f)%field%begdim3
            enddim3  = tape(t)%hlist(f)%field%enddim3
            allocate(tape(t)%hlist(f)%hbuf(begdim1:enddim1,begdim2:enddim2,begdim3:enddim3))
            tape(t)%hlist(f)%hbuf = 0._r8
            if(tape(t)%hlist(f)%field%flag_xyfill .or. (avgflag_pertape(t) .eq. 'L')) then
               allocate (tape(t)%hlist(f)%nacs(enddim1-begdim1+1,begdim3:enddim3))
            else
               allocate (tape(t)%hlist(f)%nacs(1,begdim3:enddim3))
            end if
            tape(t)%hlist(f)%nacs(:,:)=0           
            tape(t)%hlist(f)%field%vector_compliment=0            
         end do
      end do
      ! Setup vector pairs for unstructured grid interpolation
      call setup_interpolation_and_define_vector_compliments()
      !  Initialize the sat following history subsystem
      call sat_hist_init()
      
      return
   end subroutine intht


   subroutine setup_interpolation_and_define_vector_compliments()
     use dycore, only : dycore_is
     use interp_mod,       only: setup_history_interpolation, var_is_vector_uvar, var_is_vector_vvar
     use sat_hist, only : is_satfile
     integer :: t, f, ii, ff

     if(dycore_is('UNSTRUCTURED')) then
        do t=1,ptapes
           do f=1,nflds(t)     
              tape(t)%hlist(f)%field%vector_compliment=0
           end do
        end do
! Do not interpolate IC history and sat hist files
        call setup_history_interpolation(ptapes-2)
        do t=1,ptapes
           if(.not. is_satfile(t) .and. .not. is_initfile(t)) then
              do f=1,nflds(t)
                 ii=var_is_vector_uvar(tape(t)%hlist(f)%field%name)
                 if(ii>0) then
                    do ff=1,nflds(t)
                       if(ii==var_is_vector_vvar(tape(t)%hlist(ff)%field%name)) then
                          tape(t)%hlist(f)%field%vector_compliment=ff
                          tape(t)%hlist(ff)%field%vector_compliment=-f
                       end if
                    end do
                 end if
              end do
           end if
        end do
     end if
   end subroutine setup_interpolation_and_define_vector_compliments





   subroutine restart_vars_setnames()

     restartvars(1)%name = 'rgnht'
     restartvars(1)%type = pio_int
     restartvars(1)%ndims = 1     
     restartvars(1)%dims(1) = 1  ! ptapes

     restartvars(2)%name = 'nhtfrq'
     restartvars(2)%type = pio_int
     restartvars(2)%ndims = 1
     restartvars(2)%dims(1) = 1  ! ptapes


     restartvars(3)%name = 'nflds'
     restartvars(3)%type = pio_int
     restartvars(3)%ndims = 1
     restartvars(3)%dims(1) = 1  ! ptapes

     restartvars(4)%name = 'nfils'
     restartvars(4)%type = pio_int
     restartvars(4)%ndims = 1
     restartvars(4)%dims(1) = 1  ! ptapes

     restartvars(5)%name = 'mfilt'
     restartvars(5)%type = pio_int
     restartvars(5)%ndims = 1
     restartvars(5)%dims(1) = 1  ! ptapes

     restartvars(6)%name = 'nfpath'
     restartvars(6)%type = pio_char
     restartvars(6)%ndims = 2
     restartvars(6)%dims(1) = 2  ! max_string_len
     restartvars(6)%dims(2) = 1  ! ptapes

     restartvars(7)%name = 'cpath'
     restartvars(7)%type = pio_char
     restartvars(7)%ndims = 2
     restartvars(7)%dims(1) = 2  ! max_string_len
     restartvars(7)%dims(2) = 1  ! ptapes

     restartvars(8)%name = 'nhfil'
     restartvars(8)%type = pio_char
     restartvars(8)%ndims = 2
     restartvars(8)%dims(1) = 2  ! max_string_len
     restartvars(8)%dims(2) = 1  ! ptapes

     restartvars(9)%name = 'ndens'
     restartvars(9)%type = pio_int
     restartvars(9)%ndims = 1
     restartvars(9)%dims(1) = 1  ! ptapes

     restartvars(10)%name = 'fincllonlat'
     restartvars(10)%type = pio_char
     restartvars(10)%ndims = 3
     restartvars(10)%dims(1) = 5  ! max_chars
     restartvars(10)%dims(2) = 4  ! pflds
     restartvars(10)%dims(3) = 1  ! ptapes

     restartvars(11)%name = 'ncprec'
     restartvars(11)%type = pio_int
     restartvars(11)%ndims = 1
     restartvars(11)%dims(1) = 1  ! ptapes

     restartvars(12)%name = 'ngroup'
     restartvars(12)%type = pio_int
     restartvars(12)%ndims = 1
     restartvars(12)%dims(1) = 1  ! ptapes

     restartvars(13)%name = 'beg_time'
     restartvars(13)%type = pio_double
     restartvars(13)%ndims = 1
     restartvars(13)%dims(1) = 1  ! ptapes

     restartvars(14)%name = 'fincl'
     restartvars(14)%type = pio_char
     restartvars(14)%ndims = 3
     restartvars(14)%dims(1) = 3  ! fieldname_lenp2
     restartvars(14)%dims(2) = 4  ! pflds
     restartvars(14)%dims(3) = 1  ! ptapes

     restartvars(15)%name = 'fexcl'
     restartvars(15)%type = pio_char
     restartvars(15)%ndims = 3
     restartvars(15)%dims(1) = 3  ! fieldname_lenp2
     restartvars(15)%dims(2) = 4  ! pflds
     restartvars(15)%dims(3) = 1  ! ptapes

     restartvars(16)%name = 'field_name'
     restartvars(16)%type = pio_char
     restartvars(16)%ndims = 3
     restartvars(16)%dims(1) = 6  ! max_fieldname_len
     restartvars(16)%dims(2) = 7  ! maxnflds
     restartvars(16)%dims(3) = 1  ! ptapes


     restartvars(17)%name = 'decomp_type'
     restartvars(17)%type = pio_int
     restartvars(17)%ndims = 2
     restartvars(17)%dims(1) = 7  ! maxnflds
     restartvars(17)%dims(2) = 1  ! ptapes


     restartvars(18)%name = 'numlev'
     restartvars(18)%type = pio_int
     restartvars(18)%ndims = 2
     restartvars(18)%dims(1) = 7  ! maxnflds
     restartvars(18)%dims(2) = 1  ! ptapes

     restartvars(19)%name = 'hrestpath'
     restartvars(19)%type = pio_char
     restartvars(19)%ndims = 2
     restartvars(19)%dims(1) = 2  ! max_string_len
     restartvars(19)%dims(2) = 1  ! ptapes


     restartvars(20)%name = 'hwrt_prec'
     restartvars(20)%type = pio_int
     restartvars(20)%ndims = 2
     restartvars(20)%dims(1) = 7  ! maxnflds
     restartvars(20)%dims(2) = 1  ! ptapes

     restartvars(21)%name = 'avgflag'
     restartvars(21)%type = pio_char
     restartvars(21)%ndims = 2
     restartvars(21)%dims(1) = 7  ! maxnflds
     restartvars(21)%dims(2) = 1  ! ptapes

     restartvars(22)%name = 'sampling_seq'
     restartvars(22)%type = pio_char
     restartvars(22)%ndims = 3
     restartvars(22)%dims(1) = 5  ! max_chars
     restartvars(22)%dims(2) = 7  ! maxnflds
     restartvars(22)%dims(3) = 1  ! ptapes

     restartvars(23)%name = 'long_name'
     restartvars(23)%type = pio_char
     restartvars(23)%ndims = 3
     restartvars(23)%dims(1) = 5  ! max_chars
     restartvars(23)%dims(2) = 7  ! maxnflds
     restartvars(23)%dims(3) = 1  ! ptapes

     restartvars(24)%name = 'units'
     restartvars(24)%type = pio_char
     restartvars(24)%ndims = 3
     restartvars(24)%dims(1) = 5  ! max_chars
     restartvars(24)%dims(2) = 7  ! maxnflds
     restartvars(24)%dims(3) = 1  ! ptapes

     restartvars(25)%name = 'xyfill'
     restartvars(25)%type = pio_int
     restartvars(25)%ndims = 2
     restartvars(25)%dims(1) = 7  ! maxnflds
     restartvars(25)%dims(2) = 1  ! ptapes     

     restartvars(26)%name = 'lcltod_start'
     restartvars(26)%type = pio_int
     restartvars(26)%ndims = 1
     restartvars(26)%dims(1) = 1  ! ptapes     

     restartvars(27)%name = 'lcltod_stop'
     restartvars(27)%type = pio_int
     restartvars(27)%ndims = 1
     restartvars(27)%dims(1) = 1  ! ptapes     

     restartvars(28)%name = 'fillvalue'
     restartvars(28)%type = pio_double
     restartvars(28)%ndims = 2
     restartvars(28)%dims(1) = 7  ! maxnflds
     restartvars(28)%dims(2) = 1  ! ptapes     

     restartvars(29)%name = 'mdims'
     restartvars(29)%type = pio_int
     restartvars(29)%ndims = 3
     restartvars(29)%dims(1) = 8  ! maxvarmdims
     restartvars(29)%dims(2) = 7  ! maxnflds
     restartvars(29)%dims(3) = 1  ! ptapes     


   end subroutine restart_vars_setnames

   subroutine restart_dims_setnames()
     restartdims(1)%name = 'ptapes'
     restartdims(1)%len  = ptapes

     restartdims(2)%name = 'max_string_len'
     restartdims(2)%len  = max_string_len

     restartdims(3)%name = 'fieldname_lenp2'
     restartdims(3)%len  = fieldname_lenp2

     restartdims(4)%name = 'pflds'
     restartdims(4)%len  = pflds

     restartdims(5)%name = 'max_chars'
     restartdims(5)%len  = max_chars

     restartdims(6)%name = 'max_fieldname_len'
     restartdims(6)%len  = max_fieldname_len

     restartdims(7)%name = 'maxnflds'
     restartdims(7)%len  = maxval(nflds)

     restartdims(8)%name = 'maxvarmdims'
     restartdims(8)%len  = maxvarmdims

   end subroutine restart_dims_setnames


   subroutine init_restart_history (File)

!--------------------------------------------------------------------------------------------------
!
! Arguments
!
      type(file_desc_t), intent(inout) :: File                      ! Pio file Handle
!
! Local 
!
      integer :: dimids(4), ndims
      integer :: ierr, i, k

      ! Don't need to write restart data if we have written the file this step
      where (hstwr(:))
         rgnht(:) = .false.
      elsewhere
         rgnht(:) = .true.
      end where

      if(maxval(nflds)>0) then
         call restart_vars_setnames()
         call restart_dims_setnames()

	 call pio_seterrorhandling(File, PIO_BCAST_ERROR)
         do i=1,restartdimcnt
	    ! it's possible that one or more of these have been defined elsewhere
	    ierr = pio_inq_dimid(File,restartdims(i)%name, restartdims(i)%dimid)
	    if(ierr/=PIO_NOERR) then
               ierr = pio_def_dim(File,restartdims(i)%name, restartdims(i)%len, &
                    restartdims(i)%dimid)
            end if
         end do
	 call pio_seterrorhandling(File, PIO_INTERNAL_ERROR)

         do i=1,restartvarcnt
            ndims= restartvars(i)%ndims
            do k=1,ndims
               dimids(k)=restartdims(restartvars(i)%dims(k))%dimid
            end do
            allocate(restartvars(i)%vdesc)
            ierr = pio_def_var(File, restartvars(i)%name, restartvars(i)%type, dimids(1:ndims), restartvars(i)%vdesc)
         end do
      end if
    end subroutine init_restart_history

    function restartvar_getdesc(name) result(vdesc)
      character(len=*), intent(in) :: name
      type(var_desc_t), pointer :: vdesc
      character(len=max_chars) :: errmsg
      integer :: i
      
      nullify(vdesc)
      do i=1,restartvarcnt
         if(name .eq. restartvars(i)%name) then
            vdesc=>restartvars(i)%vdesc
            exit
         end if
      end do
      if(.not.associated(vdesc)) then
         errmsg = 'Could not find restart variable '//name
         call endrun(errmsg)
      end if
    end function restartvar_getdesc


!#######################################################################

    subroutine write_restart_history ( File, & 
         yr_spec, mon_spec, day_spec, sec_spec )

      implicit none
      !--------------------------------------------------------------------------------------------------
      !
      ! Arguments
      !
      type(file_desc_t), intent(inout) :: file         ! PIO restart file pointer
      integer, intent(in), optional :: yr_spec         ! Simulation year
      integer, intent(in), optional :: mon_spec        ! Simulation month
      integer, intent(in), optional :: day_spec        ! Simulation day
      integer, intent(in), optional :: sec_spec        ! Seconds into current simulation day
      !
      ! Local workspace
      !
      integer :: ierr, t, f
      integer :: rgnht_int(ptapes), start(2), count(2), startc(3)
      type(var_desc_t), pointer :: vdesc

      type(var_desc_t), pointer ::  field_name_desc, decomp_type_desc, numlev_desc, &
           hrestpath_desc, avgflag_desc, sseq_desc, longname_desc, units_desc, hwrt_prec_desc, &
           xyfill_desc, mdims_desc, fillval_desc

      integer, allocatable :: xyfill(:,:), allmdims(:,:,:)
      integer :: maxnflds


      maxnflds = maxval(nflds)
      allocate(xyfill(maxnflds, ptapes))
      xyfill = 0
      !
      !-----------------------------------------------------------------------
      ! Write the history restart data if necessary
      !-----------------------------------------------------------------------

      rgnht_int(:) = 0

      if(.not.allocated(restarthistory_tape)) allocate(restarthistory_tape(ptapes))

      do t=1,ptapes
         ! No need to write history IC restart because it is always instantaneous
         if (is_initfile(file_index=t)) rgnht(t) = .false.
         ! No need to write restart data for empty files
         if (nflds(t) == 0) rgnht(t) = .false.
         if(rgnht(t)) then
            rgnht_int(t) = 1
            restarthistory_tape(t)%hlist => history_tape(t)%hlist

            if(associated(history_tape(t)%column)) then
               restarthistory_tape(t)%column => history_tape(t)%column
            end if
            if(associated(history_tape(t)%column_st)) then
               restarthistory_tape(t)%column_st => history_tape(t)%column_st
            end if
         end if
      end do

      if(maxval(nflds)<=0) return
      
      call wshist(rgnht)

      vdesc => restartvar_getdesc('fincl')
      ierr= pio_put_var(File, vdesc, fincl(:,1:ptapes))

      vdesc => restartvar_getdesc('fincllonlat')
      ierr= pio_put_var(File, vdesc, fincllonlat(:,1:ptapes))

      vdesc => restartvar_getdesc('fexcl')
      ierr= pio_put_var(File, vdesc, fexcl(:,1:ptapes))

      vdesc => restartvar_getdesc('rgnht')
      ierr= pio_put_var(File, vdesc, rgnht_int(1:ptapes))

      vdesc => restartvar_getdesc('nhtfrq')
      ierr= pio_put_var(File, vdesc, nhtfrq(1:ptapes))

      vdesc => restartvar_getdesc('nflds')
      ierr= pio_put_var(File, vdesc, nflds(1:ptapes))

      vdesc => restartvar_getdesc('nfils')
      ierr= pio_put_var(File, vdesc, nfils(1:ptapes))

      vdesc => restartvar_getdesc('mfilt')
      ierr= pio_put_var(File, vdesc, mfilt(1:ptapes))

      vdesc => restartvar_getdesc('nfpath')
      ierr= pio_put_var(File, vdesc, nfpath(1:ptapes))

      vdesc => restartvar_getdesc('cpath')
      ierr= pio_put_var(File, vdesc,  cpath(1:ptapes))

      vdesc => restartvar_getdesc('nhfil')
      ierr= pio_put_var(File, vdesc, nhfil(1:ptapes))

      vdesc => restartvar_getdesc('ndens')
      ierr= pio_put_var(File, vdesc, ndens(1:ptapes))
      vdesc => restartvar_getdesc('ncprec')
      ierr= pio_put_var(File, vdesc, ncprec(1:ptapes))
      vdesc => restartvar_getdesc('ngroup')
      ierr= pio_put_var(File, vdesc, ngroup(1:ptapes))
      vdesc => restartvar_getdesc('beg_time')
      ierr= pio_put_var(File, vdesc, beg_time(1:ptapes))

      vdesc => restartvar_getdesc('hrestpath')
      ierr = pio_put_var(File, vdesc, hrestpath(1:ptapes))

      vdesc => restartvar_getdesc('lcltod_start')
      ierr = pio_put_var(File, vdesc, lcltod_start(1:ptapes))

      vdesc => restartvar_getdesc('lcltod_stop')
      ierr = pio_put_var(File, vdesc, lcltod_stop(1:ptapes))

      field_name_desc => restartvar_getdesc('field_name')
      decomp_type_desc => restartvar_getdesc('decomp_type')
      numlev_desc => restartvar_getdesc('numlev')
      hwrt_prec_desc => restartvar_getdesc('hwrt_prec')

      sseq_desc => restartvar_getdesc('sampling_seq')
      longname_desc => restartvar_getdesc('long_name')
      units_desc => restartvar_getdesc('units')
      avgflag_desc => restartvar_getdesc('avgflag')
      xyfill_desc => restartvar_getdesc('xyfill')

      mdims_desc => restartvar_getdesc('mdims')
      fillval_desc => restartvar_getdesc('fillvalue')

      tape=>history_tape

      allocate(allmdims(maxvarmdims,maxval(nflds),ptapes))
      allmdims=-1

      startc(1)=1
      do t = 1,ptapes
         start(2)=t
	 startc(3)=t
         do f=1,nflds(t)
            start(1)=f
	    startc(2)=f
            ierr = pio_put_var(File, field_name_desc,startc,tape(t)%hlist(f)%field%name)
            ierr = pio_put_var(File, decomp_type_desc,start,tape(t)%hlist(f)%field%decomp_type)
            ierr = pio_put_var(File, numlev_desc,start,tape(t)%hlist(f)%field%numlev)

            ierr = pio_put_var(File, hwrt_prec_desc,start,tape(t)%hlist(f)%hwrt_prec)
            ierr = pio_put_var(File, sseq_desc,startc,tape(t)%hlist(f)%field%sampling_seq)
            ierr = pio_put_var(File, longname_desc,startc,tape(t)%hlist(f)%field%long_name)
            ierr = pio_put_var(File, units_desc,startc,tape(t)%hlist(f)%field%units)
            ierr = pio_put_var(File, avgflag_desc,start, tape(t)%hlist(f)%avgflag)

            ierr = pio_put_var(File, fillval_desc,start, tape(t)%hlist(f)%field%fillvalue)
            if(associated(tape(t)%hlist(f)%field%mdims)) then
               allmdims(1:size(tape(t)%hlist(f)%field%mdims),f,t) = tape(t)%hlist(f)%field%mdims
            end if
            if(tape(t)%hlist(f)%field%flag_xyfill)             xyfill(f,t) = 1
         end do
      end do
      ierr = pio_put_var(File, xyfill_desc, xyfill)
      ierr = pio_put_var(File, mdims_desc, allmdims)

      deallocate(xyfill, allmdims)
      return

    end subroutine write_restart_history


!#######################################################################

   subroutine read_restart_history (File)
     use pio, only : pio_inquire_variable, pio_inq_dimid, pio_inquire, pio_inq_varid, &
          pio_inquire_dimension
     
      use ppgrid,        only: begchunk, endchunk
      use phys_grid,     only: get_ncols_p
      use rgrid,         only: nlon
      use dycore,        only: dycore_is
      use cam_pio_utils, only: cam_pio_openfile, phys_decomp, dyn_decomp, dyn_stagger_decomp, get_decomp
      use ioFileMod,     only: getfil
      use sat_hist,     only: sat_hist_define, sat_hist_init
    
      use shr_sys_mod,   only: shr_sys_getenv
#if (defined SPMD)
      use spmd_utils,    only : iam, mpicom, mpichar
#endif
!
!-----------------------------------------------------------------------
!
! Arguments
!
      type(file_desc_t), intent(inout) :: File            ! unit number
!
! Local workspace
!
      integer t, f                     ! tape, field indices
      integer c                        ! chunk or lat index
      integer lenc                     ! length of useful character data
      integer numlev                   ! number of vertical levels (dimension and loop)
      integer begdim2                  ! on-node vert start index
      integer enddim2                  ! on-node vert end index
      integer ioerr                    ! error code from read()
      integer begdim1                  ! on-node dim1 start index
      integer enddim1                  ! on-node dim1 end index
      integer begdim3                  ! on-node chunk or lat start index
      integer enddim3                  ! on-node chunk or lat end index
      integer ncol                     ! number of active columns per chunk
      integer lenarr                   ! global size of array to be read
      integer flag_xyfill_int


      integer rgnht_int(ptapes)
      integer :: ierr

      character(len=max_string_len)  :: locfn       ! Local filename
      character(len=max_fieldname_len), allocatable :: tmpname(:,:)
      integer, allocatable :: decomp(:,:), tmpnumlev(:,:)
      integer, pointer :: nacs(:,:)    ! accumulation counter
      character(len=max_fieldname_len) :: fname_tmp ! local copy of field name
      character(len=max_chars), pointer :: tmpstr

      real(r8) :: sum_r8
      real(r4) :: sum_r4

      type (dim_index_3d) :: dimind, dimind_xyz, dimind_xzy    ! 3-D dimension index
      integer :: ii, i, j, k, ptapes_dimid, vsize
      integer :: beglatxy, endlatxy, beglat, endlat, beglonxy, endlonxy, beglon, endlon

      type(var_desc_t) :: vdesc, longname_desc, units_desc, avgflag_desc, sseq_desc, fillval_desc
      type(io_desc_t), pointer :: iodesc
      real(r8), allocatable :: tmpfldr8(:)
      real(r4), allocatable :: tmpfldr4(:)
      integer, allocatable :: tmpfldi(:), tmpprec(:,:)
      character(len=1), allocatable:: tmpavgflag(:, :)
      integer, allocatable :: xyfill(:,:), allmdims(:,:,:)
      integer :: nacsdimcnt, nacsval
      integer :: maxnflds, dimid
      type(master_entry), pointer :: listentry


      integer levdim             ! level dimension id
      character(len=16) :: name
      integer :: ndims, dimids(8)
      integer :: tmpdims(8), dimcnt
      integer :: value
      integer :: mtapes, mdimcnt
!
! Get users logname and machine hostname
!
      if ( masterproc )then
         logname = ' '
         call shr_sys_getenv ('LOGNAME',logname,ierr)
         host = ' '
         call shr_sys_getenv ('HOST',host,ierr)
      end if
#ifdef SPMD
      ! PIO requires netcdf attributes have consistant values on all tasks
      call mpibcast(logname, len(logname), mpichar, 0, mpicom)
      call mpibcast(host, len(host), mpichar, 0, mpicom)
#endif


      beglonxy = get_dyn_grid_parm('beglonxy')
      endlonxy = get_dyn_grid_parm('endlonxy')
      beglatxy = get_dyn_grid_parm('beglatxy')
      endlatxy = get_dyn_grid_parm('endlatxy')
      beglat = get_dyn_grid_parm('beglat')
      endlat = get_dyn_grid_parm('endlat')
      beglon = get_dyn_grid_parm('beglon')
      endlon = get_dyn_grid_parm('endlon')


      call pio_seterrorhandling(File, PIO_BCAST_ERROR)

      ierr = pio_inq_dimid(File, 'ptapes', ptapes_dimid)
      if(ierr/= PIO_NOERR) then
         if(masterproc) write(iulog,*) 'Not reading history info from restart file', ierr
         return   ! no history info in restart file
      end if
      call pio_seterrorhandling(File, PIO_INTERNAL_ERROR)

      ierr = pio_inq_dimlen(File, ptapes_dimid, mtapes)

      ierr = pio_inq_dimid(File, 'maxnflds', dimid)
      ierr = pio_inq_dimlen(File, dimid, maxnflds)

      ierr = pio_inq_dimid(File, 'maxvarmdims', dimid)
      ierr = pio_inq_dimlen(File, dimid, maxvarmdims)
      

      ierr = pio_inq_varid(File, 'rgnht', vdesc)
      ierr = pio_get_var(File, vdesc, rgnht_int(1:mtapes))      

      ierr = pio_inq_varid(File, 'nhtfrq', vdesc)
      ierr = pio_get_var(File, vdesc, nhtfrq(1:mtapes))

      ierr = pio_inq_varid(File, 'nflds', vdesc)
      ierr = pio_get_var(File, vdesc, nflds(1:mtapes))
      ierr = pio_inq_varid(File, 'nfils', vdesc)
      ierr = pio_get_var(File, vdesc, nfils(1:mtapes))
      ierr = pio_inq_varid(File, 'mfilt', vdesc)
      ierr = pio_get_var(File, vdesc, mfilt(1:mtapes))

      ierr = pio_inq_varid(File, 'nfpath', vdesc)
      ierr = pio_get_var(File, vdesc, nfpath(1:mtapes))
      ierr = pio_inq_varid(File, 'cpath', vdesc)
      ierr = pio_get_var(File, vdesc, cpath(1:mtapes))
      ierr = pio_inq_varid(File, 'nhfil', vdesc)
      ierr = pio_get_var(File, vdesc, nhfil(1:mtapes))
      ierr = pio_inq_varid(File, 'hrestpath', vdesc)
      ierr = pio_get_var(File, vdesc, hrestpath(1:mtapes))


      ierr = pio_inq_varid(File, 'ndens', vdesc)
      ierr = pio_get_var(File, vdesc, ndens(1:mtapes))
      ierr = pio_inq_varid(File, 'ncprec', vdesc)
      ierr = pio_get_var(File, vdesc, ncprec(1:mtapes))
      ierr = pio_inq_varid(File, 'ngroup', vdesc)
      ierr = pio_get_var(File, vdesc, ngroup(1:mtapes))
      ierr = pio_inq_varid(File, 'beg_time', vdesc)
      ierr = pio_get_var(File, vdesc, beg_time(1:mtapes))


      ierr = pio_inq_varid(File, 'fincl', vdesc)
      ierr = pio_get_var(File, vdesc, fincl(:,1:mtapes))

      ierr = pio_inq_varid(File, 'fincllonlat', vdesc)
      ierr = pio_get_var(File, vdesc, fincllonlat(:,1:mtapes))

      ierr = pio_inq_varid(File, 'fexcl', vdesc)
      ierr = pio_get_var(File, vdesc, fexcl(:,1:mtapes))

      ierr = pio_inq_varid(File, 'lcltod_start', vdesc)
      ierr = pio_get_var(File, vdesc, lcltod_start(1:mtapes))

      ierr = pio_inq_varid(File, 'lcltod_stop', vdesc)
      ierr = pio_get_var(File, vdesc, lcltod_stop(1:mtapes))




      allocate(tmpname(maxnflds, mtapes), decomp(maxnflds, mtapes), tmpnumlev(maxnflds,mtapes))
      ierr = pio_inq_varid(File, 'field_name', vdesc)
      ierr = pio_get_var(File, vdesc, tmpname)

      ierr = pio_inq_varid(File, 'decomp_type', vdesc)
      ierr = pio_get_var(File, vdesc, decomp)
      ierr = pio_inq_varid(File, 'numlev', vdesc)
      ierr = pio_get_var(File, vdesc, tmpnumlev)

      allocate(tmpprec(maxnflds,mtapes))
      ierr = pio_inq_varid(File, 'hwrt_prec',vdesc)
      ierr = pio_get_var(File, vdesc, tmpprec(:,:))

      allocate(xyfill(maxnflds,mtapes))
      ierr = pio_inq_varid(File, 'xyfill', vdesc)
      ierr = pio_get_var(File, vdesc, xyfill)

      allocate(allmdims(maxvarmdims,maxnflds,mtapes))
      ierr = pio_inq_varid(File, 'mdims', vdesc)
      ierr = pio_get_var(File, vdesc, allmdims)


      ierr = pio_inq_varid(File, 'avgflag', avgflag_desc)      

      ierr = pio_inq_varid(File, 'long_name', longname_desc)      
      ierr = pio_inq_varid(File, 'units', units_desc)      
      ierr = pio_inq_varid(File, 'sampling_seq', sseq_desc)      

      ierr = pio_inq_varid(File, 'fillvalue', fillval_desc)      


      rgnht(:)=.false.

      allocate(history_tape(mtapes))

      tape => history_tape

      do t=1,mtapes
         nullify(tape(t)%column)
         nullify(tape(t)%column_st)

         if(rgnht_int(t)==1) rgnht(t)=.true.
         

         call strip_null(nfpath(t))
         call strip_null(cpath(t))
         call strip_null(hrestpath(t))
         allocate(tape(t)%hlist(nflds(t)))
         
         do f=1,nflds(t)
            nullify(tape(t)%hlist(f)%field%mdims)
            ierr = pio_get_var(File,fillval_desc, (/f,t/), tape(t)%hlist(f)%field%fillvalue)
            ierr = pio_get_var(File,avgflag_desc, (/f,t/), tape(t)%hlist(f)%avgflag)
            ierr = pio_get_var(File,longname_desc, (/1,f,t/), tape(t)%hlist(f)%field%long_name)
            ierr = pio_get_var(File,units_desc, (/1,f,t/), tape(t)%hlist(f)%field%units)
            tape(t)%hlist(f)%field%sampling_seq(1:max_chars) = ' '
            ierr = pio_get_var(File,sseq_desc, (/1,f,t/), tape(t)%hlist(f)%field%sampling_seq)
            call strip_null(tape(t)%hlist(f)%field%sampling_seq)
            if(xyfill(f,t) ==1) then
               tape(t)%hlist(f)%field%flag_xyfill=.true.
            else
               tape(t)%hlist(f)%field%flag_xyfill=.false.
            end if
            call strip_null(tmpname(f,t))
            tape(t)%hlist(f)%field%name = tmpname(f,t)
            tape(t)%hlist(f)%field%decomp_type = decomp(f,t)
            tape(t)%hlist(f)%field%numlev = tmpnumlev(f,t)
            tape(t)%hlist(f)%hwrt_prec = tmpprec(f,t)

            mdimcnt = count(allmdims(:,f,t)>0)
            if(mdimcnt>0) then
               allocate(tape(t)%hlist(f)%field%mdims(mdimcnt))
               tape(t)%hlist(f)%field%mdims = allmdims(1:mdimcnt,f,t)
            end if

         end do
      end do
      deallocate(tmpname, tmpnumlev, tmpprec, decomp, xyfill)

      do t=1,ptapes
         do f=1,nflds(t)
            nullify(tape(t)%hlist(f)%field%colperdim3)
            numlev   = tape(t)%hlist(f)%field%numlev
            select case (tape(t)%hlist(f)%field%decomp_type)
            case (phys_decomp)
               tape(t)%hlist(f)%field%begdim2= 1
               tape(t)%hlist(f)%field%enddim2= numlev
               tape(t)%hlist(f)%field%begdim3 = begchunk
               tape(t)%hlist(f)%field%enddim3 = endchunk
               allocate (tape(t)%hlist(f)%field%colperdim3(begchunk:endchunk))
               do c=begchunk,endchunk
                  ncol = get_ncols_p(c)
                  tape(t)%hlist(f)%field%colperdim3(c) = ncol
               end do
               tape(t)%hlist(f)%field%begdim1  = 1
               tape(t)%hlist(f)%field%enddim1  = pcols
            case (dyn_decomp)               
               tape(t)%hlist(f)%field%begdim1  = beglonxy
               tape(t)%hlist(f)%field%enddim1  = endlonxy
               tape(t)%hlist(f)%field%begdim3  = beglatxy
               tape(t)%hlist(f)%field%enddim3  = endlatxy
               tape(t)%hlist(f)%field%begdim2  = 1
               tape(t)%hlist(f)%field%enddim2  = numlev
               
               if (.not. dycore_is('LR') )then                  
                  allocate (tape(t)%hlist(f)%field%colperdim3(beglat:endlat))
                  do c=beglat,endlat
                     tape(t)%hlist(f)%field%colperdim3(c) = nlon(c)
                  end do
               end if
            case (dyn_stagger_decomp)
                  tape(t)%hlist(f)%field%begdim1 = beglonxy
                  tape(t)%hlist(f)%field%enddim1 = endlonxy
                  tape(t)%hlist(f)%field%begdim2 = 1
                  tape(t)%hlist(f)%field%enddim2 = numlev
                  tape(t)%hlist(f)%field%begdim3 = beglatxy
                  tape(t)%hlist(f)%field%enddim3 = endlatxy
            case default
               write(iulog,*)'READ_RESTART_HISTORY: bad decomp_type=',tape(t)%hlist(f)%field%decomp_type
               call endrun ()
            end select

            begdim1  = tape(t)%hlist(f)%field%begdim1
            enddim1  = tape(t)%hlist(f)%field%enddim1
            begdim2  = tape(t)%hlist(f)%field%begdim2
            enddim2  = tape(t)%hlist(f)%field%enddim2
            begdim3  = tape(t)%hlist(f)%field%begdim3
            enddim3  = tape(t)%hlist(f)%field%enddim3
            allocate(tape(t)%hlist(f)%hbuf(begdim1:enddim1,begdim2:enddim2,begdim3:enddim3))

            nullify(tape(t)%hlist(f)%varid)
            nullify(tape(t)%hlist(f)%nacs)
            if(tape(t)%hlist(f)%field%flag_xyfill .or. (avgflag_pertape(t)=='L')) then
               allocate(tape(t)%hlist(f)%nacs(enddim1-begdim1+1,begdim3:enddim3))
            else
               allocate(tape(t)%hlist(f)%nacs(1,begdim3:enddim3))
            end if
            ! initialize all buffers to zero - this will be overwritten later by the
            ! data in the history restart file if it exists.
            call h_zero(f,t)


         end do
      end do
!
!-----------------------------------------------------------------------
! Read history restart files
!-----------------------------------------------------------------------
!
! Loop over the total number of history files declared and
! read the pathname for any history restart files
! that are present (if any). Test to see if the run is a restart run
! AND if any history buffer regen files exist (rgnht=.T.). Note, rgnht 
! is preset to false, reset to true in routine WSDS if hbuf restart files
! are written and saved in the master restart file. Each history buffer
! restart file is then obtained.
! Note: some f90 compilers (e.g. SGI) complain about I/O of 
! derived types which have pointer components, so explicitly read each one.
! 
      do t=1,mtapes
         if (rgnht(t)) then
!
! Open history restart file
!
            call getfil (hrestpath(t), locfn)
            call cam_pio_openfile(tape(t)%File, locfn, 0)
!
! Read history restart file
!
            !     Read the mdims on file - they are the last dims defined and 'lev' is always the first.
            ierr=pio_inq_dimid (tape(t)%File, 'lev', levdim)
            ierr = pio_inquire(tape(t)%File, ndims)
            
            do i=levdim,ndims
               ierr = pio_inquire_dimension(tape(t)%File, i, name, value)
               call register_hist_mdim(name,value)
            enddo

            do f=1,nflds(t)  
               begdim1    =  tape(t)%hlist(f)%field%begdim1
               enddim1    =  tape(t)%hlist(f)%field%enddim1
               begdim2    =  tape(t)%hlist(f)%field%begdim2
               enddim2    =  tape(t)%hlist(f)%field%enddim2
               begdim3    =  tape(t)%hlist(f)%field%begdim3
               enddim3    =  tape(t)%hlist(f)%field%enddim3
               numlev     =  tape(t)%hlist(f)%field%numlev

               fname_tmp = strip_suffix(tape(t)%hlist(f)%field%name)
               if(masterproc) write(iulog, *) 'Reading history variable ',fname_tmp
               ierr = pio_inq_varid(tape(t)%File, fname_tmp, vdesc)


               if(.not. associated(tape(t)%hlist(f)%field%mdims)) then
                  ierr = pio_inquire_variable(tape(t)%File, vdesc, ndims=ndims,dimids=dimids)
                  dimcnt=0
                  do i=1,ndims
                     if(dimids(i)>=levdim) then
                        dimcnt=dimcnt+1
                        tmpdims(dimcnt)=dimids(i)-levdim+1
                     end if
                  end do
                  if(dimcnt>0) then
                     allocate(tape(t)%hlist(f)%field%mdims(dimcnt))
                     tape(t)%hlist(f)%field%mdims(:) = tmpdims(1:dimcnt)
                     if(dimcnt > maxvarmdims) maxvarmdims=dimcnt
                  end if
               end if
               call get_decomp(iodesc, tape(t)%hlist(f)%field, pio_double)
               call pio_setframe(vdesc, int(1,kind=PIO_OFFSET))
               call pio_read_darray(tape(t)%File, vdesc, iodesc, tape(t)%hlist(f)%hbuf, ierr)

               ierr = pio_inq_varid(tape(t)%File, trim(fname_tmp)//'_nacs', vdesc)
	       ierr = pio_inq_varndims(tape(t)%File, vdesc, nacsdimcnt)

               if(nacsdimcnt>0) then
                  call get_decomp(iodesc, tape(t)%hlist(f)%field, pio_int, numlev_in=1)
                  allocate(tape(t)%hlist(f)%nacs(begdim1:enddim1,begdim3:enddim3))

                  nacs       => tape(t)%hlist(f)%nacs(:,:)
                  call pio_read_darray(tape(t)%File, vdesc, iodesc, nacs, ierr)
               else
                  allocate(tape(t)%hlist(f)%nacs(1,begdim3:enddim3))
                  ierr = pio_get_var(tape(t)%File, vdesc, nacsval)
                  tape(t)%hlist(f)%nacs(1,:)= nacsval
               end if



            end do
!          
! Done reading this history restart file
!
            call pio_closefile (tape(t)%File)

        end if  ! rgnht(t)
        call column_init(t, nflds(t))
      end do     ! end of do mtapes loop

      


!
! If the history files are partially complete (contain less than
! mfilt(t) time samples, then get the files and open them.)
!
! NOTE:  No need to perform this operation for IC history files or empty files
!

      do t=1,mtapes
         if (is_initfile(file_index=t)) then
!
! Initialize filename specifier for IC file
!
            hfilename_spec(t) = '%c.cam2.i.%y-%m-%d-%s.nc'
            nfils(t) = 0
         else if (nflds(t) == 0) then
            nfils(t) = 0
         else
            if (nfils(t) > 0) then
               call getfil (cpath(t), locfn)
               call cam_pio_openfile(tape(t)%File, locfn, PIO_WRITE)
               call h_inquire (t)
               if(is_satfile(t)) then
                  !  Initialize the sat following history subsystem
                  call sat_hist_init()
                  call sat_hist_define(tape(t)%File)
               end if
            end if
!
! If the history file is full, close the current unit
!
            if (nfils(t) >= mfilt(t)) then
               if (masterproc) then
                  write(iulog,*)'READ_RESTART_HISTORY: nf_close(',t,')=',nhfil(t), mfilt(t)
               end if
               do f=1,nflds(t)
                  deallocate(tape(t)%hlist(f)%varid)
                  nullify(tape(t)%hlist(f)%varid)
               end do
               call pio_closefile(tape(t)%File)
               nfils(t) = 0
            end if
         end if
      end do

      ! Setup vector pairs for unstructured grid interpolation
      call setup_interpolation_and_define_vector_compliments()

      if(mtapes/=ptapes .and. masterproc) then
         write(iulog,*) ' WARNING: Restart file ptapes setting ',mtapes,' not equal to model setting ',ptapes
      end if


!
! set flag indicating h-tape contents are now defined (needed by addfld)
!
      htapes_defined = .true.      


      return
   end subroutine read_restart_history

!#######################################################################

   character(len=max_string_len) function get_hfilepath( tape )
!
!----------------------------------------------------------------------- 
! 
! Purpose: Return full filepath of history file for given tape number
! This allows public read access to the filenames without making
! the filenames public data.
!
!----------------------------------------------------------------------- 
!
  integer, intent(in) :: tape  ! Tape number

  get_hfilepath = cpath( tape )
  end function get_hfilepath

!#######################################################################

   character(len=max_string_len) function get_hist_restart_filepath( tape )
!
!----------------------------------------------------------------------- 
! 
! Purpose: Return full filepath of restart file for given tape number
! This allows public read access to the filenames without making
! the filenames public data.
!
!----------------------------------------------------------------------- 
!
  integer, intent(in) :: tape  ! Tape number

  get_hist_restart_filepath = hrestpath( tape )
  end function get_hist_restart_filepath

!#######################################################################

  integer function get_ptapes( )
!
!----------------------------------------------------------------------- 
! 
! Purpose: Return the number of tapes being used.
! This allows public read access to the number of tapes without making
! ptapes public data.
!
!----------------------------------------------------------------------- 
!
  get_ptapes = ptapes
  end function get_ptapes

  recursive function get_entry_by_name(listentry, name) result(entry)
    type(master_entry),  pointer :: listentry
    character(len=*), intent(in) :: name ! variable name
    type(master_entry), pointer :: entry
    
    if(associated(listentry)) then
       if(listentry%field%name .eq. name) then
          entry => listentry
       else
          entry=>get_entry_by_name(listentry%next_entry, name)
       end if
    else
       nullify(entry)
    end if
  end function get_entry_by_name


!#######################################################################

   subroutine fldlst ()

     use dycore, only : dycore_is
!
!----------------------------------------------------------------------- 
! 
! Purpose: Define the contents of each history file based on namelist input for initial or branch
! run, and restart data if a restart run.
!          
! Method: Use arrays fincl and fexcl to modify default history tape contents.
!         Then sort the result alphanumerically for later use by OUTFLD to
!         allow an n log n search time.
!
!---------------------------Local variables-----------------------------
!
      integer t, f                   ! tape, field indices
      integer ff                     ! index into include, exclude and fprec list
      character(len=fieldname_len) :: name ! field name portion of fincl (i.e. no avgflag separator)
      character(len=max_fieldname_len) :: mastername ! name from masterlist field
      character(len=max_chars) :: tmpcolumn_name,latlonname,latlonnamep1 ! tmp char fields
      character(len=1) :: avgflag    ! averaging flag
      character(len=1) :: prec_wrt   ! history buffer write precision flag
      character(len=6) :: prec_str
      character(len=32) :: fldname

      type (hentry) :: tmp           ! temporary used for swapping
      type (column_info) :: tmpcolumn ! temporary used for swapping

      type(master_entry), pointer :: listentry

      
      tape => history_tape
!
! First ensure contents of fincl, fexcl, and fwrtpr are all valid names
!
      do t=1,ptapes
         f = 1
         do while (f < pflds .and. fincl(f,t) /= ' ')
            name = getname (fincl(f,t))
            mastername=''
            listentry => get_entry_by_name(masterlinkedlist, name)
            if(associated(listentry)) mastername = listentry%field%name
            if (name /= mastername) then
               write(iulog,*)'FLDLST: ', trim(name), ' in fincl(', f, ') not found'
               call endrun
            end if
            f = f + 1
         end do

         f = 1
         do while (f < pflds .and. fexcl(f,t) /= ' ')
            mastername=''
            listentry => get_entry_by_name(masterlinkedlist, fexcl(f,t))
            if(associated(listentry)) mastername = listentry%field%name

            if (fexcl(f,t) /= mastername) then
               write(iulog,*)'FLDLST: ', fexcl(f,t), ' in fexcl(', f, ') not found'
               call endrun
            end if
            f = f + 1
         end do

         f = 1
         do while (f < pflds .and. fwrtpr(f,t) /= ' ')
            name = getname (fwrtpr(f,t))
            mastername=''
            listentry => get_entry_by_name(masterlinkedlist, name)
            if(associated(listentry)) mastername = listentry%field%name
            if (name /= mastername) then
               write(iulog,*)'FLDLST: ', trim(name), ' in fwrtpr(', f, ') not found'
               call endrun
            end if
            do ff=1,f-1                 ! If duplicate entry is found, stop
               if (trim(name) == trim(getname(fwrtpr(ff,t)))) then
                  write(iulog,*)'FLDLST: Duplicate field ', name, ' in fwrtpr'
                  call endrun
               end if
            end do
            f = f + 1
         end do
      end do

      nflds(:) = 0
      ! IC history file is to be created, set properties
      if(is_initfile()) then
         hfilename_spec(ptapes) = '%c.cam2.i.%y-%m-%d-%s.nc'

         ncprec(ptapes) = pio_double
         ndens (ptapes) = 1
         mfilt (ptapes) = 1
      end if




      do t=1,ptapes
!
! Add the field to the tape if specified via namelist (FINCL[1-ptapes]), or if
! it is on by default and was not excluded via namelist (FEXCL[1-ptapes]).
! Also set history buffer accumulation and output precision values according
! to the values specified via namelist (FWRTPR[1-ptapes])
! or, if not on the list, to the default values given by ndens(t).
!
         listentry => masterlinkedlist
         do while(associated(listentry))
            mastername = listentry%field%name
            call list_index (fincl(1,t), mastername, ff)

            if (ff > 0) then
               nflds(t) = nflds(t)+1
            else if ((.not. empty_htapes) .or. (is_initfile(file_index=t))) then
               call list_index (fexcl(1,t), mastername, ff)
               if (ff == 0 .and. listentry%actflag(t)) then
                  nflds(t) = nflds(t)+1
               end if
            end if
            listentry=>listentry%next_entry
         end do
      enddo
!
! Determine total number of active history tapes
!
      if (masterproc) then
         do t=1,ptapes
            if (nflds(t)  ==  0) then
               write(iulog,*)'FLDLST: Tape ',t,' is empty'
            end if
         end do
      endif
      allocate(history_tape(ptapes))
      tape=>history_tape


      do t=1,ptapes
         nullify(tape(t)%hlist)
         nullify(tape(t)%column)
         nullify(tape(t)%column_st)
! Now we have a field count and can allocate
         if(nflds(t)> 0) allocate(tape(t)%hlist(nflds(t)))

         do ff=1,nflds(t)
            nullify(tape(t)%hlist(ff)%hbuf)
            nullify(tape(t)%hlist(ff)%nacs)
            nullify(tape(t)%hlist(ff)%varid)
         end do


         nflds(t) = 0 ! recount to support array based method
         listentry => masterlinkedlist
         do while(associated(listentry))
            mastername = listentry%field%name

            call list_index (fwrtpr(1,t), mastername, ff)
            if (ff > 0) then
               prec_wrt = getflag(fwrtpr(ff,t))
            else
               prec_wrt = ' '
            end if

            call list_index (fincl(1,t), mastername, ff)

            if (ff > 0) then
               avgflag = getflag (fincl(ff,t))
               call inifld (t, listentry, avgflag,  prec_wrt)
            else if ((.not. empty_htapes) .or. (is_initfile(file_index=t))) then
               call list_index (fexcl(1,t), mastername, ff)
               if (ff == 0 .and. listentry%actflag(t)) then
                  call inifld (t, listentry, ' ', prec_wrt)
               else
                  listentry%actflag(t) = .false.
               end if
            else
               listentry%actflag(t) = .false.
            end if
            listentry=>listentry%next_entry

         end do
!
! If column output is specified make sure there are some fields defined
! for that tape
!
         if (nflds(t) .eq. 0 .and. fincllonlat(1,t) .ne. ' ') then
            write(iulog,*) 'FLDLST: Column output is specified for tape ',t,' but no fields defined for that tape.'
            call endrun()
         else
            call column_init(t, nflds(t))
         end if
!
! Specification of tape contents now complete.  Sort each list of active 
! entries for efficiency in OUTFLD.  Simple bubble sort.
!
         do f=nflds(t)-1,1,-1
            do ff=1,f

               if (tape(t)%hlist(ff)%field%name > tape(t)%hlist(ff+1)%field%name) then
            
                  tmp = tape(t)%hlist(ff)
                  tape(t)%hlist(ff  ) = tape(t)%hlist(ff+1)
                  tape(t)%hlist(ff+1) = tmp

               else if (tape(t)%hlist(ff  )%field%name == tape(t)%hlist(ff+1)%field%name) then
                  
                  write(iulog,*)'FLDLST: Duplicate field ', tape(t)%hlist(ff  )%field%name
                  write(iulog,*)'t,ff,name=',t,ff,tape(t)%hlist(ff  )%field%name
                  call endrun
            
               end if

            end do
         end do
!
! Bubble sort columns check for duplicates and rebuild field_column_names in newly sorted order
!
         if (ngroup(t) .gt. 0) then
            do i=ngroup(t)-1,1,-1
               do ff=1,i
                  latlonname=trim(tape(t)%column(ff)%lon_name) // "_" // trim(tape(t)%column(ff)%lat_name)
                  latlonnamep1=trim(tape(t)%column(ff+1)%lon_name) // "_" // trim(tape(t)%column(ff+1)%lat_name)
                  if (trim(latlonname) > trim(latlonnamep1)) then
                     tmpcolumn = tape(t)%column(ff)
                     tape(t)%column(ff) = tape(t)%column(ff+1)
                     tape(t)%column(ff+1) = tmpcolumn
                     if(dycore_is('LR')) then
                        tmpcolumn = tape(t)%column_st(ff)
                        tape(t)%column_st(ff) = tape(t)%column_st(ff+1)
                        tape(t)%column_st(ff+1) = tmpcolumn
                     end if
                  else if (trim(latlonname) == trim(latlonnamep1)) then
                     write(iulog,*)'FLDLST: Duplicate column entry for tape ',t,'duplicate column name is ',trim(latlonname)
                     call endrun
                  end if
               end do
            end do
            do f=1,nflds(t)
               do i=1,ngroup(t)
                  if (ngroup(t) .gt. 0) then
                     tape(t)%hlist(f)%field_column_name(i) = trim(tape(t)%hlist(f)%field%name) // "_" // &
                          trim(tape(t)%column(i)%lon_name) // "_" // trim(tape(t)%column(i)%lat_name)
                  else 
                     tape(t)%hlist(f)%field_column_name(1) = ' '
                  end if
               end do
            end do
         end if

         if (masterproc) then
            if (nflds(t) > 0) then
               write(iulog,*) ' '
               write(iulog,*)'FLDLST: History file ', t, ' contains ', nflds(t), ' fields'

               if (is_initfile(file_index=t)) then
                  write(iulog,*) ' Write frequency:                 ',inithist,' (INITIAL CONDITIONS)'
               else if (nflds(t) > 0) then
                  if (nhtfrq(t) == 0 .and. t == 1) then
                     write(iulog,*) ' Write frequency:                  MONTHLY'
                  else
                     write(iulog,*) ' Write frequency:                 ',nhtfrq(t)
                  end if
               end if

               write(iulog,*) ' Filename specifier:              ', trim(hfilename_spec(t))

               prec_str = 'double'
               if (ndens(t) == 2) prec_str = 'single'
               write(iulog,*) ' Output precision:                ', prec_str
               write(iulog,*) ' Number of time samples per file: ', mfilt(t)

               write(iulog,*)' Included fields are:'
            end if
   
            do f=1,nflds(t)
               if (ngroup(t) .gt. 0) then
                  if (f .eq. 1) write(iulog,*) '   Fields on this tape will be output as column data (FIELD_LON_LAT)'
                  do i = 1, ngroup(t)
                     ff = (f-1)*ngroup(t) + i
                     fldname = tape(t)%hlist(f)%field_column_name(i)
                     write(iulog,9000) ff, fldname, tape(t)%hlist(f)%field%units, tape(t)%hlist(f)%field%numlev, &
                             tape(t)%hlist(f)%avgflag, tape(t)%hlist(f)%field%long_name

                  end do
               else
                  fldname = tape(t)%hlist(f)%field%name
                  write(iulog,9000) f, fldname, tape(t)%hlist(f)%field%units, tape(t)%hlist(f)%field%numlev, &
                             tape(t)%hlist(f)%avgflag, tape(t)%hlist(f)%field%long_name

               end if

            end do

         end if

      end do    ! do t=1,ptapes

9000  format(i5, 1x, a32, 1x, a16, 1x, i4, 1x, a1, 2x, 256a)
!
! Packing density, ndens: With netcdf, only 1 (nf_double) and 2 (pio_real)
! are allowed
!
      do t=1,ptapes
         if (ndens(t) == 1) then
            ncprec(t) = pio_double
         else if (ndens(t) == 2) then
            ncprec(t) = pio_real
         else
            call endrun ('FLDLST: ndens must be 1 or 2')
         end if

      end do
!
! set flag indicating h-tape contents are now defined (needed by addfld)
!
      htapes_defined = .true.      
!
!  Now that masterlinkedlist is defined, construct primary and secondary hashing
!  tables.
!
      call bld_outfld_hash_tbls()
      call bld_htapefld_indices()

      return
   end subroutine fldlst

!#######################################################################
   subroutine inifld (t, listentry, avgflag, prec_wrt)
!
!----------------------------------------------------------------------- 
! 
! Purpose: Add a field to the active list for a history tape
! 
! Method: Copy the data from the master field list to the active list for the tape
!         Also: define mapping arrays from (col,chunk) -> (lon,lat)
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------


!
! Arguments
!
      integer, intent(in) :: t            ! history tape index

      type(master_entry), pointer :: listentry

      character*1, intent(in) :: avgflag  ! averaging flag
      character*1, intent(in) :: prec_wrt ! history output precision flag
!
! Local workspace
!
      integer :: i                  ! index
      integer :: n                  ! field index on defined tape


!
! Ensure that it is not to late to add a field to the history tape
!
      if (htapes_defined) then
         call endrun ('INIFLD: Attempt to add field '//listentry%field%name//' after history files set')
      end if


      nflds(t) = nflds(t) + 1
      n = nflds(t)
!
! Copy field info.
!
      if(n > size(tape(t)%hlist)) then
         write(iulog,*) 'tape field miscount error ', n, size(tape(t)%hlist)
         call endrun()
      end if

      tape(t)%hlist(n)%field = listentry%field

      select case (prec_wrt)
      case (' ')
         if (ndens(t) == 1) then
            tape(t)%hlist(n)%hwrt_prec = 8
         else
            tape(t)%hlist(n)%hwrt_prec = 4
         end if
      case ('4')
         tape(t)%hlist(n)%hwrt_prec = 4
         if (masterproc) then
            write(iulog,*) 'INIFLD: Output data type for ', tape(t)%hlist(n)%field%name, &
                       ' is real*4'
         end if
      case ('8')
         tape(t)%hlist(n)%hwrt_prec = 8
         if (masterproc) then
            write(iulog,*) 'INIFLD: Output data type for ', tape(t)%hlist(n)%field%name, &
                       ' is real*8'
         end if
      case default
         call endrun ('INIFLD: unknown prec_wrt='//prec_wrt)
      end select
!
! Override the default averaging (masterlist) averaging flag if non-blank
!
      if (avgflag == ' ') then
         tape(t)%hlist(n)%avgflag = listentry%avgflag(t)
         tape(t)%hlist(n)%time_op = listentry%time_op(t)
      else
         tape(t)%hlist(n)%avgflag = avgflag
         select case (avgflag)
         case ('A')
            tape(t)%hlist(n)%time_op = 'mean'
         case ('B')
            tape(t)%hlist(n)%time_op = 'mean00z'
         case ('I')
            tape(t)%hlist(n)%time_op = ' '
         case ('X')
            tape(t)%hlist(n)%time_op = 'maximum'
         case ('M')
            tape(t)%hlist(n)%time_op = 'minimum'
         case('L')
            tape(t)%hlist(n)%time_op = LT_DESC
         case default
            call endrun ('INIFLD: unknown avgflag='//avgflag)
         end select
      end if

#ifdef HDEBUG
      write(iulog,*)'HDEBUG: ',__LINE__,' field ', tape(t)%hlist(n)%field%name, ' added as ', 'field number ', n,' on tape ', t
      write(iulog,*)'units=',tape(t)%hlist(n)%field%units
      write(iulog,*)'numlev=',tape(t)%hlist(n)%field%numlev
      write(iulog,*)'avgflag=',tape(t)%hlist(n)%avgflag
      write(iulog,*)'time_op=',tape(t)%hlist(n)%time_op
      write(iulog,*)'hwrt_prec=',tape(t)%hlist(n)%hwrt_prec
#endif

      return
   end subroutine inifld


   subroutine column_init(t, nflds)
     use dyn_grid, only: get_dyn_grid_parm_real2d, get_dyn_grid_parm_real1d
     use dycore, only : dycore_is
     use spmd_utils, only : mpi_sum, mpi_integer, iam, mpicom
     use phys_grid, only : get_ncols_p, get_lat_all_p, get_lon_all_p
     use ppgrid,        only: begchunk, endchunk

     integer, intent(in) :: t, nflds
     integer :: ff, i

     real(r8), allocatable :: lontmp   (:)    ! Temp array holding longitude values 
     real(r8), allocatable :: lontmp_st(:)   ! Temp array holding longitude values 
     character(len=max_chars) :: lonlatname(pflds), lonname(pflds), latname(pflds) ! variable name
     integer :: lonind   (pflds,2) !   beginning and ending longitude range
     integer :: latind   (pflds,2) !   beginning and ending latitude range

     integer :: lonind_st(pflds,2) !   beginning and ending staggered longitude range
     integer :: latind_st(pflds,2) !   beginning and ending staggeredlatitude range
     real(r8), pointer :: londeg(:,:), londeg_st(:,:), latdeg(:), latdeg_st(:)
     character(len=max_chars) :: tmpname(pflds) ! variable name (will not be used)

     integer :: group              !   number of group
     integer :: splon, plat, plon
     integer, pointer :: hmask(:,:)
     integer :: lats(pcols), lons(pcols)
     integer :: cnt, cnttot, ierr, lchnk, ncol, j
     integer :: beglatxy,endlatxy,beglonxy,endlonxy
     

     plon = get_dyn_grid_parm('plon')
     splon = get_dyn_grid_parm('splon')
     plat = get_dyn_grid_parm('plat')

     beglonxy = get_dyn_grid_parm('beglonxy')
     beglatxy = get_dyn_grid_parm('beglatxy')
     endlonxy = get_dyn_grid_parm('endlonxy')
     endlatxy = get_dyn_grid_parm('endlatxy')


     londeg => get_dyn_grid_parm_real2d('londeg')
     latdeg => get_dyn_grid_parm_real1d('latdeg')
     if(dycore_is('LR')) then
        londeg_st => get_dyn_grid_parm_real2d('londeg_st')
        latdeg_st => get_dyn_grid_parm_real1d('latdeg_st')
     end if
     !
     ! Setup column information if this field will be written as group
     ! First verify the column information in the namelist
     !

     ff = 1

     allocate(lontmp(plon))
     if(dycore_is('LR')) then
        allocate(lontmp_st(splon))
     end if
     do while (fincllonlat(ff,t) /= ' ')

        lonlatname(ff) = trim(fincllonlat(ff,t))
        call getlatind(lonlatname(ff),latind(ff,1),latind(ff,2),latname(ff),latdeg     ,plat)
	
        do i = 1,plon
           lontmp(i) = londeg(i,1)
           if (lontmp(i) .lt. 0._r8) lontmp(i) = lontmp(i) + 360._r8
        end do
        call getlonind(lonlatname(ff),lonind(ff,1),lonind(ff,2),lonname(ff),lontmp,plon)
#ifdef HDEBUG
        write(iulog,*)'HDEBUG: ',__LINE__,' closest lon lat range for group ',lonlatname(ff),' is ',lonind(ff,:),latind(ff,:)
        write(iulog,*) londeg(lonind(ff,1),1),londeg(lonind(ff,2),1),latdeg(latind(ff,1)),latdeg(latind(ff,2))
#endif
        if(dycore_is('LR')) then
           call getlatind(lonlatname(ff),latind_st(ff,1),latind_st(ff,2),tmpname(ff),latdeg_st     ,plat-1)
           do i = 1,splon
              lontmp_st(i) = londeg_st(i,1)
              if (lontmp_st(i) .lt. 0._r8) lontmp_st(i) = lontmp_st(i) + 360._r8
           end do
           call getlonind(lonlatname(ff),lonind_st(ff,1),lonind_st(ff,2),tmpname(ff),lontmp_st,splon )
#ifdef HDEBUG
           write(iulog,*)'HDEBUG: ',__LINE__,' closest staggered lon lat range for group ',lonlatname(ff),' is ',lonind_st(ff,:),latind_st(ff,:)
           write(iulog,*) londeg_st(lonind_st(ff,1),1),londeg_st(lonind_st(ff,2),1),latdeg_st(latind_st(ff,1)),latdeg_st(latind_st(ff,2))
#endif
        end if
        ff = ff + 1
     end do
     deallocate(lontmp)

     if(dycore_is('LR')) then      
        deallocate(lontmp_st)
     endif

     group = ff-1
     ngroup(t)=0

     if (group .gt. 0) then 
        if(collect_column_output(t)) then
           allocate(tape(t)%column(1))
           tape(t)%column(1)%num_lats = 0
           tape(t)%column(1)%num_lons = 0

           allocate(hmask(plon, plat))
         
           hmask = 0
           cnt = 0
           do j=1,plat
              do i=1,plon
                 do ff=1,group
                    if(hmask(i,j)==0 .and. i>=lonind(ff,1) .and. i <= lonind(ff,2) .and. &
                       j>=latind(ff,1) .and. j <= latind(ff,2)) then
                       cnt = cnt+1
                       hmask(i,j)=cnt
                    end if
                 end do                 
              end do
           end do

           tape(t)%column(1)%column_cnt = cnt
           allocate(tape(t)%column(1)%hmask(pcols,begchunk:endchunk))
           tape(t)%column(1)%hmask=0
           allocate(tape(t)%column(1)%hmask_dyn(endlonxy-beglonxy+1,endlatxy-beglatxy+1))
           tape(t)%column(1)%hmask_dyn=hmask(beglonxy:endlonxy,beglatxy:endlatxy)

           do lchnk=begchunk,endchunk
              ncol = get_ncols_p(lchnk)
              call get_lat_all_p(lchnk, pcols, lats)
              call get_lon_all_p(lchnk, pcols, lons)
              do i=1,ncol
                 tape(t)%column(1)%hmask(i,lchnk)=hmask(lons(i),lats(i))
              end do
           end do

           deallocate(hmask)


        else
           ngroup(t)=group
           allocate (tape(t)%column   (group))
           if(dycore_is('LR')) then
              allocate (tape(t)%column_st(group))
           endif
           do i=1,nflds
              allocate (tape(t)%hlist(i)%field_column_name(group))
           end do
           do ff = 1, group
              nullify(tape(t)%column(ff)%hmask)
              tape(t)%column(ff)%lat_name=trim(latname(ff))
              tape(t)%column(ff)%lon_name=trim(lonname(ff))
              tape(t)%column(ff)%columnlon(:)=lonind(ff,:)
              tape(t)%column(ff)%columnlat(:)=latind(ff,:)
              tape(t)%column(ff)%num_lats = &
                   tape(t)%column(ff)%columnlat(2)-tape(t)%column(ff)%columnlat(1)+1
              tape(t)%column(ff)%num_lons = &
                   tape(t)%column(ff)%columnlon(2)-tape(t)%column(ff)%columnlon(1)+1
              do i=1,nflds
                 tape(t)%hlist(i)%field_column_name(ff) = &
                      trim(tape(t)%hlist(i)%field%name) // "_" // trim(lonname(ff)) // "_" // trim(latname(ff))
              end do
              if(dycore_is('LR')) then
                 tape(t)%column_st(ff)%lat_name=trim(latname(ff)) // "_st"
                 tape(t)%column_st(ff)%lon_name=trim(lonname(ff)) // "_st"
                 tape(t)%column_st(ff)%columnlon(:)=lonind_st(ff,:)
                 tape(t)%column_st(ff)%columnlat(:)=latind_st(ff,:)
                 tape(t)%column_st(ff)%num_lats = &
                      tape(t)%column_st(ff)%columnlat(2)-tape(t)%column_st(ff)%columnlat(1)+1
                 tape(t)%column_st(ff)%num_lons = &
                      tape(t)%column_st(ff)%columnlon(2)-tape(t)%column_st(ff)%columnlon(1)+1
              endif
           
           end do
        end if

     else

        allocate (tape(t)%column(1))
        nullify(tape(t)%column(1)%hmask)
        do i=1,nflds
           allocate (tape(t)%hlist(i)%field_column_name(1))
           tape(t)%hlist(i)%field_column_name(1) = ' '
        end do
        tape(t)%column(1)%lat_name=' '
        tape(t)%column(1)%lon_name=' '
        tape(t)%column(1)%columnlon(:)=0
        tape(t)%column(1)%columnlat(:)=0
        tape(t)%column(1)%num_lats=0
        tape(t)%column(1)%num_lons=0
     end if


   end subroutine column_init








!#######################################################################

   subroutine strip_null(str)
     character(len=*), intent(inout) :: str
     do i=1,len(str)
        if(ichar(str(i:i))==0) str(i:i)=' '
     end do
   end subroutine strip_null

   character(len=max_fieldname_len) function strip_suffix (name)
!
!---------------------------------------------------------- 
! 
! Purpose:  Strip "&IC" suffix from fieldnames if it exists
!          
!----------------------------------------------------------
!
! Arguments
!
      character(len=*), intent(in) :: name
!
! Local workspace
!
      integer :: n
!
!-----------------------------------------------------------------------
!
      strip_suffix = ' '

      do n = 1,fieldname_len
         strip_suffix(n:n) = name(n:n)
         if(name(n+1:n+1         ) == ' '                       ) return
         if(name(n+1:n+fieldname_suffix_len) == fieldname_suffix) return
      end do

      strip_suffix(fieldname_len+1:max_fieldname_len) = name(fieldname_len+1:max_fieldname_len)

      return

   end function strip_suffix

!#######################################################################

   character(len=fieldname_len) function getname (inname)
!
!----------------------------------------------------------------------- 
! 
! Purpose: retrieve name portion of inname
!          
! Method:  If an averaging flag separater character is present (":") in inname, 
!          lop it off
! 
!-------------------------------------------------------------------------------
!
! Arguments
!
      character(len=*), intent(in) :: inname
!
! Local workspace
!
      integer :: length
      integer :: i
   
      length = len (inname)
      
      if (length < fieldname_len .or. length > fieldname_lenp2) then
         write(iulog,*) 'GETNAME: bad length=',length
         call endrun
      end if
   
      getname = ' '
      do i=1,fieldname_len
         if (inname(i:i) == ':') exit
         getname(i:i) = inname(i:i)
      end do
      
      return
   end function getname

!#######################################################################

   subroutine getlatind(inname,beglatind,endlatind,latname,latdeg,mlat)
!
!----------------------------------------------------------------------- 
! 
! Purpose: retrieve closest model lat index given north south latitude string
!          
! Author: John Truesdale
! 
!-------------------------------------------------------------------------------
!
! Arguments
!
      character(len=*), intent(in)  :: inname
      character(len=*), intent(out) :: latname
      integer, intent(out) :: beglatind,endlatind
      integer, intent(in) :: mlat
      real(r8), intent(in) :: latdeg(mlat)
!
! Local workspace
!
      character(len=max_chars) str1,str2,tmpstr(4)
      integer :: i,j,marker,latind(2),tmplen(4)
      real(r8) :: latdegree(2)

!
!-------------------------------------------------------------------------------
!
      str1=' '
      str2=' '
      tmpstr(:)= ' '

!
! make sure _ separator is present
!      
      if (scan(inname,'_').eq.0) then
         write(iulog,*)'GETLATIND: Improperly formatted column string.  Missing underscore character (xxxE_yyyS) ', &
                   inname,scan(inname,'_')
         call endrun
      end if
!
! split inname string into lat lon strings
!      
      marker=scan(inname,'_')
      str1=inname(:marker-1)
      str2=inname(marker+1:)

!
! split ranges of lats( or lons) into seperate substrings. Substrings 1,2 will contain lats(lons) from portion of string before underscore
! if a single column and not a range is specified substrings 1 and 2 will be the same
! and substrings 3,4 will contain lats(lons) from portion of main string after underscore character
! if a single column and not a range is specified substrings 3 and 4 will be the same
!
      if (scan(str1,':') .ne. 0) then
         marker=scan(str1,':')
         tmpstr(1)=str1(:marker-1)
         tmpstr(2)=str1(marker+1:)
      else
         tmpstr(1)=str1
         tmpstr(2)=str1
      end if
      if (scan(str2,':') .ne. 0) then
         marker=scan(str2,':')
         tmpstr(3)=str2(:marker-1)
         tmpstr(4)=str2(marker+1:)
      else
         tmpstr(3)=str2
         tmpstr(4)=str2
      end if

! check format of substrings - (number followed by single character north/south east/west designation)

      do i = 1,4
         tmplen(i)=len_trim(tmpstr(i))
         if (verify(tmpstr(i),"0123456789.").ne.tmplen(i) .or. verify(tmpstr(i)(tmplen(i):tmplen(i)),"eEwWnNsS") .ne. 0) then
            write(iulog,*)'GETLATIND (2): Improperly formatted column string. ',&
                 inname,verify(tmpstr(i),"0123456789."),tmplen(i),&
                 verify(tmpstr(i)(tmplen(i):tmplen(i)),"eEwWnNsS")
            call endrun
         end if
      end do

! find latitude substrings and put into temporary work space tmpstr(1:2)

      if (verify(tmpstr(1)(tmplen(1):tmplen(1)),"nNsSs").eq.0  .and. verify(tmpstr(2)(tmplen(2):tmplen(2)),"nNsSs").eq.0 ) then

      else if (verify(tmpstr(3)(tmplen(3):tmplen(3)),"nNsSs").eq.0 .and. verify(tmpstr(3)(tmplen(3):tmplen(3)),"nNsSs").eq.0) then
         tmpstr(1) = tmpstr(3)
         tmplen(1)=tmplen(3)
         tmpstr(2) = tmpstr(4)
         tmplen(2)=tmplen(4)
      else
         call endrun ('GETLATIND (3): Improperly formatted column string. '//inname)
      end if
!
! convert lat substrings to real and if in southern hemisphere make negative
!
      do i = 1,2
         read(tmpstr(i)(1:tmplen(i)-1),*) latdegree(i)
         if (verify(tmpstr(i)(tmplen(i):tmplen(i)),'sS').eq.0) then
            latdegree(i)=(-1._r8)*latdegree(i)
         end if
!
! Make sure specified latitudes is in bounds
!
         if (latdegree(i) .lt. -90._r8 .or. latdegree(i) .gt. 90._r8) then
            write(iulog,*)'GETLATIND: latitude for column namelist is out of range (-90 .. 90) value=',latdegree(i)
            call endrun
         endif
!
! Find closest lat index for each substring
!
         latind(i)=mlat
         do j = 1, mlat-1
            if ( abs(latdeg(j)-latdegree(i)) .lt.  &
                 abs(latdeg(j+1)-latdegree(i))) then
               latind(i)=j
               exit
            endif
         enddo
      end do
!
! output begining and ending latitude indicies.   If just a column is specified then beginning latitude index will be the same as the
! ending latitude index  

!      
      if (latind(1) .le. latind(2) ) then
         beglatind =latind(1)
         endlatind =latind(2)
      else 
         beglatind = latind(2)
         endlatind = latind(1)
      end if
      if (beglatind .eq. endlatind) then
         latname = 'LAT_'//trim(tmpstr(1))
      else
         if (latind(1) .le. latind(2) ) then
            latname = 'LAT_'//trim(tmpstr(1)) // "_to_" // trim(tmpstr(2))
         else 
            latname = 'LAT_'//trim(tmpstr(2)) // "_to_" // trim(tmpstr(1))
         end if
      end if
      ! one last sanity check
      if(scan(latname,' ') > 0) then
         latname(scan(latname,' '):) = ' '
      end if

      return
   end subroutine getlatind
!#######################################################################

   subroutine getlonind (inname,beglonind,endlonind,lonname,londeg,mlon)
!
!----------------------------------------------------------------------- 
! 
! Purpose: retrieve closes model lat index given north south latitude string
!          
! Method:  If an averaging flag separater character is present (":") in inname, 
!          lop it off
! 
!-------------------------------------------------------------------------------
!
! Arguments
!
      character(len=*) inname,lonname
      integer :: beglonind,endlonind
      integer :: mlon
      real(r8) :: londeg(mlon)
!
! Local workspace
!
      character(len=max_chars) str1,str2,tmpstr(4)
      integer :: i,j,marker,lonind(2),tmplen(4)
      real(r8) :: londegree(2)
      real(r8) :: min, tmp
!
!-------------------------------------------------------------------------------
!
      str1=' '
      str2=' '
      tmpstr(:)= ' '
!
! make sure _ separator is present
!      
      if (scan(inname,'_').eq.0) then
         write(iulog,*)'GETLONIND: Improperly formatted column string.  Missing underscore character (xxxE_yyyS) ',&
              inname,scan(inname,'_')
         call endrun
      end if
!
! split string in to lat lon strings
!      
      marker=scan(inname,'_')
      str1=inname(:marker-1)
      str2=inname(marker+1:)

!
! split ranges of lats( or lons) into seperate substrings. Substrings 1,2 will contain lats(lons) from portion of string before underscore
! if a single column and not a range is specified substrings 1 and 2 will be the same
! and substrings 3,4 will contain lats(lons) from portion of main string after underscore character
! if a single column and not a range is specified substrings 3 and 4 will be the same
!
      if (scan(str1,':') .ne. 0) then
         marker=scan(str1,':')
         tmpstr(1)=str1(:marker-1)
         tmpstr(2)=str1(marker+1:)
      else
         tmpstr(1)=str1
         tmpstr(2)=str1
      end if

      if (scan(str2,':') .ne. 0) then
         marker=scan(str2,':')
         tmpstr(3)=str2(:marker-1)
         tmpstr(4)=str2(marker+1:)
      else
         tmpstr(3)=str2
         tmpstr(4)=str2
      end if

! check format of substrings - (number followed by single character north/south east/west designation)

      do i = 1,4
         tmplen(i)=len_trim(tmpstr(i))
         if (verify(tmpstr(i),"0123456789.").ne.tmplen(i) .or. verify(tmpstr(i)(tmplen(i):tmplen(i)),"eEwWnNsS") .ne. 0) then
            write(iulog,*)'GETLONIND (2): Improperly formatted column string. ', &
                 inname,verify(tmpstr(i),"0123456789."),tmplen(i), &
                 verify(tmpstr(i)(tmplen(i):tmplen(i)),"eEwWnNsS")
            call endrun
         end if
      end do

! find longitude substrings and put into temporary work space tmpstr(1:2)

      if (verify(tmpstr(1)(tmplen(1):tmplen(1)),"eEwW").eq.0  .and. verify(tmpstr(2)(tmplen(2):tmplen(2)),"eEwW").eq.0 ) then

      else if (verify(tmpstr(3)(tmplen(3):tmplen(3)),"eEwW").eq.0 .and. verify(tmpstr(3)(tmplen(3):tmplen(3)),"eEwW").eq.0) then
         tmpstr(1) = tmpstr(3)
         tmplen(1)=tmplen(3)
         tmpstr(2) = tmpstr(4)
         tmplen(2)=tmplen(4)
      else
         call endrun ('GETLONIND (3): Improperly formatted column string. '//inname)
      end if
!
! convert lon substrings to real and make sure its degrees east
!
      do i = 1,2
         read(tmpstr(i)(1:tmplen(i)-1),*) londegree(i)
         if (verify(tmpstr(i)(tmplen(i):tmplen(i)),'wW').eq.0) then
            londegree(i) = 360._r8 - londegree(i)
         end if
!
! Make sure specified longitudes are in bounds
!
         if (londegree(i) .lt. 0._r8 .or. londegree(i) .gt. 360._r8) then
            write(iulog,*)'GETLONIND: longitude for column namelist is out of range (0 .. 360) value=',londegree(i)
            call endrun
         endif
!
! Find closest lon index for each substring.  If just a column is specified then beginning longitude index will be the same as the
! ending longitude index  
!
         min  = 1.e+36_r8
         do j = 1, mlon
            tmp = abs(londeg(j)-londegree(i))
            if ( tmp .le. min) then
               min       = tmp
               lonind(i) = j
            endif
         end do
      end do
!
! output begining and ending longitude indicies.   If just a column is specified then beginning longitude index will be the same as the
! ending longitude index  
!      
      if (lonind(1) .le. lonind(2) ) then
         beglonind =lonind(1)
         endlonind =lonind(2)
      else 
         beglonind = lonind(2)
         endlonind = lonind(1)
      end if


      if (beglonind .eq. endlonind) then
         lonname = 'LON_'//trim(tmpstr(1))
      else
         if (lonind(1) .le. lonind(2) ) then
            lonname = 'LON_'//trim(tmpstr(1)) // "_to_" // trim(tmpstr(2))
         else 
            lonname = 'LON_'//trim(tmpstr(2)) // "_to_" // trim(tmpstr(1))
         end if
      end if


      return
    end subroutine getlonind

!#######################################################################

   character(len=1) function getflag (inname)
!
!----------------------------------------------------------------------- 
! 
! Purpose: retrieve flag portion of inname
!          
! Method:  If an averaging flag separater character is present (":") in inname, 
!          return the character after it as the flag
! 
!-------------------------------------------------------------------------------
!
! Arguments
!
      character(len=*), intent(in) :: inname   ! character string
!
! Local workspace
!
      integer :: length         ! length of inname
      integer :: i              ! loop index

      length = len (inname)

      if (length /= fieldname_lenp2) then
         write(iulog,*) 'GETFLAG: bad length=',length
         call endrun
      end if

      getflag = ' '
      do i=1,fieldname_lenp2-1
         if (inname(i:i) == ':') then
            getflag = inname(i+1:i+1)
            exit
         end if
      end do

      return
   end function getflag

!#######################################################################

   subroutine list_index (list, name, index)
!
! Input arguments
!
      character(len=*) , intent(in) :: list(pflds) ! input list of names, possibly ":" delimited
      character(len=max_fieldname_len), intent(in) :: name ! name to be searched for
!
! Output arguments
!
      integer, intent(out) :: index               ! index of "name" in "list"
!
! Local workspace
!
      character(len=fieldname_len) :: listname    ! input name with ":" stripped off.
      integer f                       ! field index

      index = 0
      do f=1,pflds
!
! Only list items
!
         listname = getname (list(f))
         if (listname == ' ') exit
         if (listname == name) then
            index = f
            exit
         end if
      end do
      
      return
   end subroutine list_index

!#######################################################################

   subroutine outfld (fname, field, idim, c)
     use cam_history_buffers, only: hbuf_accum_inst, hbuf_accum_add, &
          hbuf_accum_add00z, hbuf_accum_max, hbuf_accum_min, &
          hbuf_accum_addlcltime, dim_index_2d
     
!
!----------------------------------------------------------------------- 
! 
! Purpose: Accumulate (or take min, max, etc. as appropriate) input field
!          into its history buffer for appropriate tapes
! 
! Method: Check 'masterlist' whether the requested field 'fname' is active
!         on one or more history tapes, and if so do the accumulation.
!         If not found, return silently.
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
!
! Arguments
!
      character(len=*), intent(in) :: fname ! Field name--should be 8 chars long

      integer, intent(in) :: idim           ! Longitude dimension of field array
      integer, intent(in) :: c              ! chunk (physics) or latitude (dynamics) index

      real(r8), intent(in) :: field(idim,*) ! Array containing field values
!
! Local variables
!
      integer :: t, f                ! tape, field indices
      integer :: fl, fu              ! upper, lower indices used in binary search thru sorted list
      integer :: begdim1             ! on-node dim1 start index
      integer :: enddim1             ! on-node dim1 end index
      integer :: begdim2             ! on-node dim2 start index
      integer :: enddim2             ! on-node dim2 end index
      integer :: begdim3
      integer :: enddim3
      integer :: endi                ! ending longitude index (reduced grid)

      character*(max_fieldname_len) :: fname_loc  ! max-char equivalent of fname
      character*1 :: avgflag         ! averaging flag
      
      real(r8),pointer :: hbuf(:,:)      ! history buffer
      integer, pointer :: nacs(:)    ! accumulation counter
      type (dim_index_2d) :: dimind  ! 2-D dimension index
      logical :: flag_xyfill         ! non-applicable xy points flagged with fillvalue
      integer :: ff                  ! masterlist index pointer
      integer :: ierr, ncol
!-----------------------------------------------------------------------

      tape => history_tape

      fname_loc = fname

      ff = get_masterlist_indx(fname_loc)

!
!  If ( ff < 0 ), the field is not defined on the masterlist. This check
!  is necessary because of coding errors calling outfld without first defining
!  the field on masterlist.
!
      if ( ff < 0 ) then
         return
      end if
!
!  Next, check to see whether this field is active on one or more history
!  tapes.
!
      if ( .not. masterlist(ff)%thisentry%act_sometape )  then
         return
      end if
!
! Note, the field may be on any or all of the history files (primary
! and auxiliary).
!
!      write(iulog,*)'fname_loc=',fname_loc
      do 40 t=1,ptapes
         if ( .not. masterlist(ff)%thisentry%actflag(t)) cycle
         f = masterlist(ff)%thisentry%htapeindx(t)
!
! Update history buffer
!
         begdim1 = tape(t)%hlist(f)%field%begdim1
         enddim1 = tape(t)%hlist(f)%field%enddim1
         begdim2 = tape(t)%hlist(f)%field%begdim2
         enddim2 = tape(t)%hlist(f)%field%enddim2
         begdim3 = tape(t)%hlist(f)%field%begdim3
         enddim3 = tape(t)%hlist(f)%field%enddim3

         avgflag = tape(t)%hlist(f)%avgflag
         flag_xyfill = tape(t)%hlist(f)%field%flag_xyfill

         nacs   => tape(t)%hlist(f)%nacs(:,c)

         
         hbuf => tape(t)%hlist(f)%hbuf(:,:,c)

         if(associated(tape(t)%hlist(f)%field%colperdim3)) then         
            endi    = tape(t)%hlist(f)%field%colperdim3(c)
         else
            endi = tape(t)%hlist(f)%field%enddim1 - tape(t)%hlist(f)%field%begdim1 + 1
         end if
         dimind = dim_index_2d (1,endi,begdim2,enddim2)

         select case (avgflag)

         case ('I') ! Instantaneous
            call hbuf_accum_inst (hbuf, field, nacs, dimind, idim, flag_xyfill,tape(t)%hlist(f)%field%fillvalue)

         case ('A') ! Time average

            call hbuf_accum_add (hbuf, field, nacs, dimind, idim, flag_xyfill,tape(t)%hlist(f)%field%fillvalue)

         case ('B') ! Time average only 00z values

            call hbuf_accum_add00z (hbuf, field, nacs, dimind, idim, flag_xyfill,tape(t)%hlist(f)%field%fillvalue)

         case ('X') ! Maximum over time

            call hbuf_accum_max (hbuf, field, nacs, dimind, idim, flag_xyfill,tape(t)%hlist(f)%field%fillvalue)

         case ('M') ! Minimum over time

           call hbuf_accum_min (hbuf, field, nacs, dimind, idim, flag_xyfill,tape(t)%hlist(f)%field%fillvalue)

        case ('L')

           call hbuf_accum_addlcltime(hbuf, field, nacs, dimind, idim, flag_xyfill, &
                tape(t)%hlist(f)%field%fillvalue, &
                c, tape(t)%hlist(f)%field%decomp_type, &
                lcltod_start(t),lcltod_stop(t))

        case default

            call endrun ('OUTFLD: invalid avgflag='//avgflag)

         end select
40    continue

      return
   end subroutine outfld

!#######################################################################

   logical function is_initfile (file_index)
!
!------------------------------------------------------------------------ 
! 
! Purpose: to determine:
!
!   a) if an IC file is active in this model run at all
!       OR,
!   b) if it is active, is the current file index referencing the IC file
!      (IC file is always at ptapes)
! 
!------------------------------------------------------------------------
!
! Arguments
!
      integer, intent(in), optional :: file_index ! index of file in question

      is_initfile = .false.

      if (present(file_index)) then
         if (inithist /= 'NONE' .and. file_index == ptapes) is_initfile = .true.
      else
         if (inithist /= 'NONE'                           ) is_initfile = .true.
      end if

      return

   end function is_initfile

!#######################################################################

   integer function strcmpf (name1, name2)
!
!----------------------------------------------------------------------- 
! 
! Purpose: Return the lexical difference between two strings
! 
! Method: Use ichar() intrinsic as we loop through the names
! 
!-----------------------------------------------------------------------
!
! Arguments
!
      character(len=max_fieldname_len), intent(in) :: name1, name2 ! strings to compare
      integer n                                     ! loop index
   
      do n=1,max_fieldname_len
         strcmpf = ichar(name1(n:n)) - ichar(name2(n:n))
         if (strcmpf /= 0) exit
      end do

      return
   end function strcmpf

!#######################################################################

   subroutine h_inquire (t)
     use pio, only : pio_inquire, pio_inq_dimid, pio_inquire_dimension, pio_inq_varid, &
          pio_inquire_variable, pio_inq_attlen
     use dycore, only : dycore_is
!
!----------------------------------------------------------------------- 
! 
! Purpose: Ensure that the proper variables are on a history file
! 
! Method: Issue the appropriate netcdf wrapper calls
! 
!-----------------------------------------------------------------------
!
! Arguments
!
      integer, intent(in) :: t   ! tape index
!
! Local workspace
!
      integer f,ff               ! field index
      integer ret                ! return value from function call


      integer :: ierr
      integer :: len
      integer :: mdimsize

!
!
! Dimension id's
!
      tape => history_tape


      
!
! Create variables for model timing and header information 
!
      if(.not. is_satfile(t)) then
         ierr=pio_inq_varid (tape(t)%File,'ndcur   ',    tape(t)%ndcurid)
         ierr=pio_inq_varid (tape(t)%File,'nscur   ',    tape(t)%nscurid)
         ierr=pio_inq_varid (tape(t)%File,'nsteph  ',    tape(t)%nstephid)

         ierr=pio_inq_varid (tape(t)%File,'time_bnds',   tape(t)%tbndid)
         ierr=pio_inq_varid (tape(t)%File,'date_written',tape(t)%date_writtenid)
         ierr=pio_inq_varid (tape(t)%File,'time_written',tape(t)%time_writtenid)
#if ( defined BFB_CAM_SCAM_IOP )
         ierr=pio_inq_varid (tape(t)%File,'tsec    ',tape(t)%tsecid)
         ierr=pio_inq_varid (tape(t)%File,'bdate   ',tape(t)%bdateid)
#endif
         if (.not. is_initfile(file_index=t) ) then
            ! Don't write the GHG/Solar forcing data to the IC file.  It is never
            ! read from that file so it's confusing to have it there.
            ierr=pio_inq_varid (tape(t)%File,'co2vmr  ',    tape(t)%co2vmrid)
            ierr=pio_inq_varid (tape(t)%File,'ch4vmr  ',    tape(t)%ch4vmrid)
            ierr=pio_inq_varid (tape(t)%File,'n2ovmr  ',    tape(t)%n2ovmrid)
            ierr=pio_inq_varid (tape(t)%File,'f11vmr  ',    tape(t)%f11vmrid)
            ierr=pio_inq_varid (tape(t)%File,'f12vmr  ',    tape(t)%f12vmrid)
            ierr=pio_inq_varid (tape(t)%File,'sol_tsi ',    tape(t)%sol_tsiid)
         end if
      end if
      ierr=pio_inq_varid (tape(t)%File,'date    ',    tape(t)%dateid)
      ierr=pio_inq_varid (tape(t)%File,'datesec ',    tape(t)%datesecid)
      ierr=pio_inq_varid (tape(t)%File,'time    ',    tape(t)%timeid)


!
! Obtain variable name from ID which was read from restart file
!
      do f=1,nflds(t)
         if(.not.associated(tape(t)%hlist(f)%varid)) then
            allocate(tape(t)%hlist(f)%varid(max(1,ngroup(t))))
         end if
!
! If this field will be put out as columns then get column names for field
!
         if (ngroup(t) .gt. 0) then
            do i=1,ngroup(t)
               ierr=pio_inq_varid (tape(t)%File, tape(t)%hlist(f)%field_column_name(i), tape(t)%hlist(f)%varid(i))
               ierr=pio_get_att(tape(t)%File, tape(t)%hlist(f)%varid(i), 'basename',tape(t)%hlist(f)%field%name)
            end do
         else   
            ierr=pio_inq_varid (tape(t)%File,tape(t)%hlist(f)%field%name, tape(t)%hlist(f)%varid(1))
         end if
         if(tape(t)%hlist(f)%field%numlev>1) then
            ierr = pio_inq_attlen(tape(t)%File,tape(t)%hlist(f)%varid(1),'mdims', mdimsize)
            allocate(tape(t)%hlist(f)%field%mdims(mdimsize))
            ierr=pio_get_att(tape(t)%File,tape(t)%hlist(f)%varid(1),'mdims', &
                 tape(t)%hlist(f)%field%mdims(1:mdimsize))
            if(mdimsize>maxvarmdims) maxvarmdims=mdimsize
         end if

      end do

      if(masterproc) then
         write(iulog,*)'H_INQUIRE: Successfully opened netcdf file '
      end if

      return
   end subroutine h_inquire

!#######################################################################

   subroutine add_default (name, tindex, flag)
!
!----------------------------------------------------------------------- 
! 
! Purpose: Add a field to the default "on" list for a given history file
! 
! Method: 
! 
!-----------------------------------------------------------------------
!
! Arguments
!
      character(len=*), intent(in) :: name  ! field name
      character(len=1), intent(in) :: flag  ! averaging flag

      integer, intent(in) :: tindex         ! history tape index
!
! Local workspace
!
      integer :: t            ! file index
      integer :: f            ! field index
      logical :: found        ! flag indicates field found in masterlist
      type(master_entry), pointer :: listentry

!
! Check validity of input arguments
!
      if (tindex > ptapes) then
         write(iulog,*)'ADD_DEFAULT: tape index=', tindex, ' is too big'
         call endrun
      end if

! Add to IC file if tindex = 0, reset to ptapes
      if (tindex == 0) then
         t = ptapes
         if ( .not. is_initfile(file_index=t) ) return
      else
         t = tindex
      end if

      if (flag /= ' ' .and. flag /= 'A' .and. flag /= 'B' .and. flag /= 'I' .and. &
          flag /= 'X' .and. flag /= 'M' .and. flag /= 'L') then

         call endrun ('ADD_DEFAULT: unknown averaging flag='//flag)
      end if
!
! Look through master list for input field name.  When found, set active
! flag for that tape to true.  Also set averaging flag if told to use other
! than default.
!
      found = .false.
      listentry => get_entry_by_name(masterlinkedlist, trim(name))
      if(.not.associated(listentry)) then
         call endrun ('ADD_DEFAULT: field='//name//' not found')
      end if
      listentry%actflag(t) = .true.
      if (flag /= ' ') then
         listentry%avgflag(t) = flag
         select case (flag)
         case ('A')
            listentry%time_op(t) = 'mean'
         case ('B')
            listentry%time_op(t) = 'mean00z'
         case ('I')
            listentry%time_op(t) = ' '
         case ('X')
            listentry%time_op(t) = 'maximum'
         case ('M')
            listentry%time_op(t) = 'minimum'
         case ('L')
            listentry%time_op(t) = LT_DESC
         case default
            call endrun ('ADD_DEFAULT: unknown avgflag='//flag)
         end select
      end if
      
      return
   end subroutine add_default

!#######################################################################

   subroutine h_override (t)
!
!----------------------------------------------------------------------- 
! 
! Purpose: Override default history tape contents for a specific tape
!
! Method: Copy the flag into the master field list
! 
!-----------------------------------------------------------------------
!
! Arguments
!
      integer, intent(in) :: t         ! history tape index
!
! Local workspace
!
      character(len=1) :: avgflg       ! lcl equiv of avgflag_pertape(t) (to address xlf90 compiler bug)

      type(master_entry), pointer :: listentry


      avgflg = avgflag_pertape(t)

      
      listentry=>masterlinkedlist
      do while(associated(listentry))
         listentry%avgflag(t) = avgflag_pertape(t)
         select case (avgflg)
         case ('A')
            listentry%time_op(t) = 'mean'
         case ('B')
            listentry%time_op(t) = 'mean00z'
         case ('I')
            listentry%time_op(t) = ' '
         case ('X')
            listentry%time_op(t) = 'maximum'
         case ('M')
            listentry%time_op(t) = 'minimum'
         case ('L')
            listentry%time_op(t) = LT_DESC
            
         case default
            call endrun ('H_OVERRIDE: unknown avgflag='//avgflag_pertape(t))
         end select
         listentry=>listentry%next_entry
      end do

   end subroutine h_override
         
!#######################################################################

   subroutine h_define (t, restart)
     !
     !----------------------------------------------------------------------- 
     ! 
     ! Purpose: Define contents of history file t
     ! 
     ! Method: Issue the required netcdf wrapper calls to define the history file contents
     ! 
     !-----------------------------------------------------------------------
     use pspect, only : ptrn, ptrk, ptrm
     use rgrid, only : nlon, wnummax

     use dyn_grid, only : get_dyn_grid_parm_real2d, get_dyn_grid_parm_real1d
     use dycore, only : dycore_is
     use time_manager, only: get_step_size, get_ref_date, timemgr_get_calendar_cf
     use filenames,    only: caseid
     use string_utils, only: to_upper
     use abortutils,   only: endrun
     use dycore,       only: dycore_is
     use physconst,     only: pi
     use spmd_utils,    only: iam, nsmps, npes, mpicom
     use cam_pio_utils, only: cam_pio_createfile, get_phys_decomp
     use interp_mod,    only: get_interp_lat, get_interp_lon, latlon_interpolation, add_interp_attributes
     use cosp_share,   only: srmid_cosp, dbzemid_cosp, scol_cosp, sza_cosp, htmisrmid_cosp, prsmid_cosp, &
          taumid_cosp, htmid_cosp, taumid_cosp_modis, prslim_cosp, taulim_cosp, taulim_cosp_modis, &
          dbzelim_cosp, srlim_cosp, htmisrlim_cosp, htlim_cosp
     use sat_hist,     only: sat_hist_define
     
     use ppgrid, only : begchunk,endchunk, pver, pverp
     use phys_grid, only : get_rlat_all_p, get_rlon_all_p     

     !-----------------------------------------------------------------------

     !
     ! Input arguments
     !
     integer, intent(in) :: t   ! tape index
     logical, intent(in) :: restart
     !
     ! Local workspace
     !
     integer :: i, j            ! longitude, latitude indices
     integer :: f               ! field index
     integer :: ff              ! varid index for fields output by column
     integer :: numlev          ! number of vertical levels (dimension and loop)
     integer :: ncreal          ! real data type for output
     integer :: dtime           ! timestep size
     integer :: ndbase = 0      ! days component of base time
     integer :: nsbase = 0      ! seconds component of base time
     integer :: nbdate          ! base date in yyyymmdd format
     integer :: bdate           ! base date in yyyymmdd format
     integer :: nbsec           ! time of day component of base date [seconds]
     integer :: yr, mon, day    ! year, month, day components of a date


     integer :: slatdim         ! staggered latitude dimension
     integer :: slondim         ! staggered longitude dimension
     type(var_desc_t) :: slatvar         ! variable id for staggered lat
     type(var_desc_t) :: slonvar         ! variable id for staggered lon
     type(var_desc_t) :: wsid            ! Staggered latitude weight ID

     real(r8), pointer:: londeg_st(:,:)    ! Staggered grid point array (lon)
     real(r8), pointer:: latdeg_st(: )   ! Staggered grid point array (lat)
     real(r8), pointer:: w_staggered(:)  ! Staggered location weights

     real(r8) :: ailev(pverp)      ! interface level values
     real(r8), pointer :: latdeg(:)    ! degrees gaussian latitudes 
     real(r8), pointer :: latdeg2d(:,:)    ! degrees latitudes for ncols
     real(r8), pointer :: londeg(:,:)    ! degrees longitude
     real(r8), pointer :: w(:)    ! grid weights
     real(r8), parameter :: radtodeg = 180.0_r8/pi
     real(r8), pointer :: harea(:)

     character(len=max_chars) str ! character temporary 
     character(len=max_chars) :: fname_tmp ! local copy of field name
     character(len=max_chars) calendar             ! Calendar type
     !
     ! netcdf variables
     !
     integer ret                ! function return value
     integer :: timdim             ! unlimited dimension id
     type(var_desc_t) :: latvar             ! latitude variable id
     type(var_desc_t) :: lonvar             ! longitude variable id
     type(var_desc_t) :: areavar             ! cell area variable id
     type(var_desc_t), pointer :: glatvar(:)     ! column latitude variable id
     type(var_desc_t), pointer :: glonvar(:)     ! column longitude variable id
     type(var_desc_t) :: ps0var             ! variable id for PS0
     integer :: chardim            ! character dimension id

     real(r8) :: alev(pver)        ! level values (pascals)
     real(r8), allocatable :: alon(:)        ! longitude values (degrees)
     integer  :: dimenchar(2)       ! character dimension ids


     integer :: nacsdims(2)        ! dimension ids for nacs (used in restart file)
     integer :: hdims, tdims
     integer londim             ! longitude dimension id
     integer latdim             ! latitude dimension id
     integer, allocatable:: grouplondim(:) ! longitude dimension id
     integer, allocatable:: grouplatdim(:) ! latitude dimension id
     integer levdim             ! level dimension id
     integer ilevdim            ! interface dimension id

     integer bnddim            ! bounds dimension id
     type(var_desc_t):: levvar             ! level variable id
     type(var_desc_t):: ilevvar            ! intfc variable id

     type(var_desc_t) :: hyaiid             ! hybrid A coef. intfc var id
     type(var_desc_t) :: hybiid             ! hybrid B coef. intfc var id
     type(var_desc_t) :: hyamid             ! hybrid A coef. level var id
     type(var_desc_t) :: hybmid             ! hybrid B coef. level var id
     type(var_desc_t) :: ntrmid             ! M truncation parameter var id
     type(var_desc_t) :: ntrnid             ! N truncation parameter var id
     type(var_desc_t) :: ntrkid             ! K truncation parameter var id



     integer :: dimindex(8)         ! dimension ids for variable declaration

     integer, pointer :: mdims(:)
     integer :: mdimsize


     type(var_desc_t), pointer :: glatvar_st(:)     ! column staggered latitude variable id
     type(var_desc_t), pointer :: glonvar_st(:)     ! column staggered longitude variable id
     integer, allocatable :: grouplondim_st(:) ! staggered longitude dimension id
     integer, allocatable :: grouplatdim_st(:) ! staggered latitude dimension id
     integer :: splon, plon

     integer ndims
     integer dim1s,dim2s        ! global size of the first and second horizontal dim.
     integer ncol
     integer ncoldim            ! horizontal dim for non rectangular grids
     integer :: ierr
     type(var_desc_t), pointer :: varid

     integer, pointer :: ldof(:), mdimids(:)
     !
     integer :: amode

     type(var_desc_t) :: cosp_prs_var      	! COSP ISCCP mean pressure variable id
     type(var_desc_t) :: cosp_prs_bnds      	! COSP ISCCP mean pressure variable bounds
     type(var_desc_t) :: cosp_tau_var      	! COSP ISCCP mean optical depth variable id
     type(var_desc_t) :: cosp_tau_bnds      	! COSP ISCCP mean optical depth varible bounds
     type(var_desc_t) :: cosp_tau_modis_var    ! COSP MODIS mean optical depth variable id
     type(var_desc_t) :: cosp_tau_modis_bnds    ! COSP MODIS mean optical depth variable bounds
     type(var_desc_t) :: cosp_ht_var      	! COSP mean ht variable id
     type(var_desc_t) :: cosp_ht_bnds      	! COSP mean ht variable bounds
     type(var_desc_t) :: cosp_dbze_var         ! COSP mean dbze variable id
     type(var_desc_t) :: cosp_dbze_bnds         ! COSP mean dbze variable bounds
     type(var_desc_t) :: cosp_sr_var   	! COSP mean sr variable id
     type(var_desc_t) :: cosp_sr_bnds   	! COSP mean sr variable bounds
     type(var_desc_t) :: cosp_scol_var         ! COSP mean scol variable id
     type(var_desc_t) :: cosp_htmisr_var   	! COSP mean htmisr variable id
     type(var_desc_t) :: cosp_htmisr_bnds   	! COSP mean htmisr variable bounds
     type(var_desc_t) :: cosp_sza_var   	! COSP mean sza variable id

     type(io_desc_t), pointer :: iodesc

     integer :: lchnk, plat


     if(restart) then
        tape => restarthistory_tape
        if(masterproc) write(iulog,*)'Opening netcdf history restart file ', trim(hrestpath(t))
     else
        tape => history_tape
        if(masterproc) write(iulog,*)'Opening netcdf history file ', trim(nhfil(t))
     end if

     amode = PIO_CLOBBER


     if(restart) then
        call cam_pio_createfile (tape(t)%File, hrestpath(t), amode)
     else
        call cam_pio_createfile (tape(t)%File, nhfil(t), amode)
     end if
     if(is_satfile(t)) then
        ierr = pio_def_dim (tape(t)%File, 'ncol', pio_unlimited, timdim)  

        ierr=pio_def_var (tape(t)%File,'lat',pio_double,(/timdim/),latvar)
        ierr=pio_put_att (tape(t)%File, latvar, 'long_name', 'latitude')
        ierr=pio_put_att (tape(t)%File, latvar, 'units', 'degrees_north')

        ierr=pio_def_var (tape(t)%File,'lon',pio_double,(/timdim/),lonvar)
        ierr=pio_put_att (tape(t)%File, lonvar,'long_name','longitude')
        ierr=pio_put_att (tape(t)%File, lonvar,'units','degrees_east')

     else
        !
        ! Setup netcdf file - create the dimensions of lat,lon,time,level
        !     
        interpolate = (latlon_interpolation(t) .and. .not. restart )
        ncoloutput = ((dycore_is('UNSTRUCTURED') .and. .not. interpolate).or. associated(tape(t)%column(1)%hmask) .and. .not. restart)

        call get_horiz_grid_dim_d(dim1s,dim2s)

        ! First define the horizontal grid dims
        
        latdim = -1
        londim = -1
        ncoldim = -1

        plat = get_dyn_grid_parm('plat')

        if(dycore_is('UNSTRUCTURED')) then
           if(interpolate) then
              dim1s=get_dyn_grid_parm('nlon')
              dim2s=get_dyn_grid_parm('nlat')

              ret = pio_def_dim (tape(t)%File, 'lat', dim2s, latdim)
              ret = pio_def_dim (tape(t)%File, 'lon', dim1s, londim)            
              
              call add_interp_attributes(tape(t)%File)

           else
              ncol = dim1s*dim2s
              ret = pio_def_dim (tape(t)%File, 'ncol', ncol, ncoldim)
           end if
           ret = pio_put_att(tape(t)%File, PIO_GLOBAL, 'np', get_dyn_grid_parm('np'))
           ret = pio_put_att(tape(t)%File, PIO_GLOBAL, 'ne', get_dyn_grid_parm('ne'))
        else if(ncoloutput) then
           ncol = tape(t)%column(1)%column_cnt
           ret =  pio_def_dim (tape(t)%File, 'ncol', ncol, ncoldim)
        else if(ngroup(t)==0 .or. restart) then
           ret = pio_def_dim (tape(t)%File, 'lat', dim2s, latdim)
           ret = pio_def_dim (tape(t)%File, 'lon', dim1s, londim)
           if(dycore_is('LR')) then
              ret = pio_def_dim (tape(t)%File, 'slat', plat-1, slatdim)
              splon = get_dyn_grid_parm('splon')
              ret = pio_def_dim (tape(t)%File, 'slon', splon, slondim)
           end if
        endif
        

        ! Define the unlimited time dim
        ret = pio_def_dim (tape(t)%File, 'time', pio_unlimited, timdim)
        ret = pio_def_dim (tape(t)%File, 'nbnd', 2, bnddim)
        ret = pio_def_dim (tape(t)%File, 'chars', 8, chardim)

        !
        ! create dimensions for groups of contiguous columns
        ! variables with for single columns are created later in this routine
        !

        if (.not. restart .and. ngroup(t).ne.0) then
           allocate(grouplatdim(ngroup(t)), grouplondim(ngroup(t)))
           allocate(glatvar(ngroup(t)), glonvar(ngroup(t)))
           do i = 1, ngroup(t)

              if(scan(trim(tape(t)%column(i)%lat_name),' ')> 0) then
                 tape(t)%column(i)%lat_name(scan(tape(t)%column(i)%lat_name, ' '):)=' '
              end if
              if(scan(trim(tape(t)%column(i)%lon_name),' ')> 0) then
                 tape(t)%column(i)%lon_name(scan(tape(t)%column(i)%lon_name, ' '):)=' '
              end if

              ret = pio_def_dim (tape(t)%File, trim(tape(t)%column(i)%lat_name), tape(t)%column(i)%num_lats, grouplatdim(i))
              ierr=pio_def_var (tape(t)%File,tape(t)%column(i)%lat_name,pio_double,grouplatdim(i:i),glatvar(i))
              ierr=pio_put_att (tape(t)%File, glatvar(i),'long_name','latitude')
              ierr=pio_put_att (tape(t)%File, glatvar(i),'units','degrees_north')

              ret = pio_def_dim (tape(t)%File, trim(tape(t)%column(i)%lon_name), tape(t)%column(i)%num_lons, grouplondim(i))
              ierr=pio_def_var (tape(t)%File,tape(t)%column(i)%lon_name,pio_double,grouplondim(i:i),glonvar(i))
              ierr=pio_put_att (tape(t)%File, glonvar(i),'long_name','longitude')
              ierr=pio_put_att (tape(t)%File, glonvar(i),'units','degrees_east')
           end do
           if(dycore_is('LR')) then
              !
              ! For staggered (FV) fields, include staggered dimension on h-file
              !
              do f=1,nflds(t)
                 fname_tmp = tape(t)%hlist(f)%field%name
                 if(fname_tmp .eq. 'US' .or. fname_tmp .eq. 'FU_S') then
                    if(.not.allocated(grouplatdim_st)) then
                       allocate(grouplatdim_st(ngroup(t)))
                       allocate(glatvar_st(ngroup(t)))

                       do i = 1, ngroup(t)

                          if(scan(trim(tape(t)%column_st(i)%lat_name),' ')> 0) then
                             tape(t)%column_st(i)%lat_name(scan(tape(t)%column_st(i)%lat_name, ' '):)=' '
                          end if
                          if(scan(trim(tape(t)%column_st(i)%lon_name),' ')> 0) then
                             tape(t)%column_st(i)%lon_name(scan(tape(t)%column_st(i)%lon_name, ' '):)=' '
                          end if

                          ret = pio_def_dim (tape(t)%File, tape(t)%column_st(i)%lat_name, &
                               tape(t)%column_st(i)%num_lats, grouplatdim_st(i))
                          ierr=pio_def_var (tape(t)%File,tape(t)%column_st(i)%lat_name,pio_double, &
                               grouplatdim_st(i:i),glatvar_st(i))
                          ierr=pio_put_att (tape(t)%File, glatvar_st(i),'long_name','staggered latitude')
                          ierr=pio_put_att (tape(t)%File, glatvar_st(i),'units','degrees_north')
                       end do
                    end if
                 else if(fname_tmp .eq. 'VS' .or. fname_tmp .eq. 'FV_S') then
                    if(.not.allocated(grouplondim_st)) then
                       allocate(grouplondim_st(ngroup(t)))
                       allocate(glonvar_st(ngroup(t)))
                       do i = 1, ngroup(t)
                          if(scan(trim(tape(t)%column_st(i)%lat_name),' ')> 0) then
                             tape(t)%column_st(i)%lat_name(scan(tape(t)%column_st(i)%lat_name, ' '):)=' '
                          end if
                          if(scan(trim(tape(t)%column_st(i)%lon_name),' ')> 0) then
                             tape(t)%column_st(i)%lon_name(scan(tape(t)%column_st(i)%lon_name, ' '):)=' '
                          end if
                          ret = pio_def_dim (tape(t)%File, tape(t)%column_st(i)%lon_name, &
                               tape(t)%column_st(i)%num_lons, grouplondim_st(i))
                          ierr=pio_def_var (tape(t)%File,tape(t)%column_st(i)%lon_name,pio_double, &
                               grouplondim_st(i:i),glonvar_st(i))
                          ierr=pio_put_att (tape(t)%File, glonvar_st(i),'long_name','staggered longitude')
                          ierr=pio_put_att (tape(t)%File, glonvar_st(i),'units','degrees_east')
                       end do
                    end if
                 end if
              end do
           endif
        end if
     endif   ! is satfile
     ! Restart requires that the def_dim call for registeredmdims be the last in the history file.



     ! Now define the registered mdims
     allocate(mdimids(registeredmdims))
     do i=1,registeredmdims
        ret = pio_def_dim (tape(t)%File, trim(hist_mdims(i)%name), hist_mdims(i)%value, mdimids(i))        
        ! define variables to label the dimensions, use the same names as the
        ! dimensions

        select case(trim(hist_mdims(i)%name))
        case ('lev')
           levdim = mdimids(i)
           ierr=pio_def_var (tape(t)%File,'lev',pio_double,(/levdim/),levvar)
           str = 'hybrid level at midpoints (1000*(A+B))'
           ierr=pio_put_att (tape(t)%File, levvar, 'long_name', trim(str))
           str = 'level'
           ierr=pio_put_att (tape(t)%File, levvar, 'units', trim(str))
           ierr=pio_put_att (tape(t)%File, levvar, 'positive', 'down')
           ierr=pio_put_att (tape(t)%File, levvar, 'standard_name', 'atmosphere_hybrid_sigma_pressure_coordinate')
           ierr=pio_put_att (tape(t)%File, levvar, 'formula_terms', 'a: hyam b: hybm p0: P0 ps: PS')

           ierr=pio_def_var (tape(t)%File,'hyam',PIO_DOUBLE,(/levdim/),hyamid)
           str = 'hybrid A coefficient at layer midpoints'
           ierr=pio_put_att (tape(t)%File, hyamid, 'long_name', trim(str))

           ierr=pio_def_var (tape(t)%File,'hybm',PIO_DOUBLE,(/levdim/),hybmid)
           str = 'hybrid B coefficient at layer midpoints'
           ierr=pio_put_att (tape(t)%File, hybmid, 'long_name', trim(str))

        case ('ilev')
           ilevdim = mdimids(i)
           ierr=pio_def_var (tape(t)%File,'ilev',pio_double,(/ilevdim/),ilevvar)
           str = 'hybrid level at interfaces (1000*(A+B))'
           ierr=pio_put_att (tape(t)%File, ilevvar, 'long_name', trim(str))
           str = 'level'
           ierr=pio_put_att (tape(t)%File, ilevvar, 'units', trim(str))
           ierr=pio_put_att (tape(t)%File, ilevvar, 'positive', 'down')
           ierr=pio_put_att (tape(t)%File, ilevvar, 'standard_name', 'atmosphere_hybrid_sigma_pressure_coordinate')
           ierr=pio_put_att (tape(t)%File, ilevvar, 'formula_terms', 'a: hyai b: hybi p0: P0 ps: PS')

           ierr=pio_def_var (tape(t)%File,'hyai',PIO_DOUBLE,(/ilevdim/),hyaiid)
           str = 'hybrid A coefficient at layer interfaces'
           ierr=pio_put_att (tape(t)%File, hyaiid, 'long_name', trim(str))

           ierr=pio_def_var (tape(t)%File,'hybi',PIO_DOUBLE,(/ilevdim/),hybiid)
           str = 'hybrid B coefficient at layer interfaces'
           ierr=pio_put_att (tape(t)%File, hybiid, 'long_name', trim(str))

	end select
        if(.not. is_satfile(t) .or. is_initfile(t)) then
           select case(trim(hist_mdims(i)%name))
           case('cosp_prs')
              ierr=pio_def_var (tape(t)%File, 'cosp_prs', PIO_DOUBLE,  mdimids(i:i), cosp_prs_var)
              str = 'COSP Mean ISCCP pressure'
              ierr=pio_put_att (tape(t)%File, cosp_prs_var, 'long_name', trim(str))
              str = 'mb'
              ierr=pio_put_att (tape(t)%File, cosp_prs_var, 'units', trim(str))
              str = 'cosp_prs_bnds'
              ierr=pio_put_att (tape(t)%File, cosp_prs_var, 'bounds', trim(str))
              ierr=pio_def_var (tape(t)%File, str, PIO_DOUBLE, (/bnddim,mdimids(i:i)/), cosp_prs_bnds)

           case('cosp_tau')
              ierr=pio_def_var (tape(t)%File, 'cosp_tau', PIO_DOUBLE,  mdimids(i:i), cosp_tau_var)
              str = 'COSP Mean ISCCP optical depth'
              ierr=pio_put_att (tape(t)%File, cosp_tau_var, 'long_name', trim(str))
              str = 'unitless'
              ierr=pio_put_att (tape(t)%File, cosp_tau_var, 'units', trim(str))
              str = 'cosp_tau_bnds'
              ierr=pio_put_att (tape(t)%File, cosp_tau_var, 'bounds', trim(str))
              ierr=pio_def_var (tape(t)%File, str, PIO_DOUBLE, (/bnddim,mdimids(i:i)/), cosp_tau_bnds)
           case('cosp_tau_modis')
              ierr=pio_def_var (tape(t)%File, 'cosp_tau_modis', PIO_DOUBLE,  mdimids(i:i), cosp_tau_modis_var)
              str = 'COSP Mean MODIS optical depth'
              ierr=pio_put_att (tape(t)%File, cosp_tau_modis_var, 'long_name', trim(str))
              str = 'unitless'
              ierr=pio_put_att (tape(t)%File, cosp_tau_modis_var, 'units', trim(str))
              str = 'cosp_tau_modis_bnds'
              ierr=pio_put_att (tape(t)%File, cosp_tau_modis_var, 'bounds', trim(str))
              ierr=pio_def_var (tape(t)%File, str, PIO_DOUBLE, (/bnddim,mdimids(i:i)/), cosp_tau_modis_bnds)
           case('cosp_ht')
              ierr=pio_def_var (tape(t)%File, 'cosp_ht', PIO_DOUBLE, mdimids(i:i), cosp_ht_var)
              str = 'COSP Mean Height for lidar and radar simulator outputs'
              ierr=pio_put_att (tape(t)%File, cosp_ht_var, 'long_name', trim(str))
              str = 'm'
              ierr=pio_put_att (tape(t)%File, cosp_ht_var, 'units', trim(str))
              str = 'cosp_ht_bnds'
              ierr=pio_put_att (tape(t)%File, cosp_ht_var, 'bounds', trim(str))
              ierr=pio_def_var (tape(t)%File, str, PIO_DOUBLE, (/bnddim,mdimids(i:i)/), cosp_ht_bnds)
           case('cosp_dbze')
              ierr=pio_def_var (tape(t)%File, 'cosp_dbze', PIO_DOUBLE,  mdimids(i:i),cosp_dbze_var)
              str = 'COSP Mean dBZe for radar simulator CFAD output'
              ierr=pio_put_att (tape(t)%File, cosp_dbze_var, 'long_name', trim(str))
              str = 'dbze'
              ierr=pio_put_att (tape(t)%File, cosp_dbze_var, 'units', trim(str))
              str = 'cosp_dbze_bnds'
              ierr=pio_put_att (tape(t)%File, cosp_dbze_var, 'bounds', trim(str))
              ierr=pio_def_var (tape(t)%File, str, PIO_DOUBLE, (/bnddim,mdimids(i:i)/), cosp_dbze_bnds)
           case('cosp_sr')
              ierr=pio_def_var (tape(t)%File, 'cosp_sr', PIO_DOUBLE,  mdimids(i:i),cosp_sr_var)
              str = 'COSP Mean Scattering Ratio for lidar simulator CFAD output'
              ierr=pio_put_att (tape(t)%File, cosp_sr_var, 'long_name', trim(str))
              str = '1'
              ierr=pio_put_att (tape(t)%File, cosp_sr_var, 'units', trim(str))
              str = 'cosp_sr_bnds'
              ierr=pio_put_att (tape(t)%File, cosp_sr_var, 'bounds', trim(str))
              ierr=pio_def_var (tape(t)%File, str, PIO_DOUBLE, (/bnddim,mdimids(i:i)/), cosp_sr_bnds)           
           case('cosp_scol')
              ierr=pio_def_var (tape(t)%File, 'cosp_scol', PIO_DOUBLE,  mdimids(i:i),cosp_scol_var)
              str = 'COSP subcolumn'
              ierr=pio_put_att (tape(t)%File, cosp_scol_var, 'long_name', trim(str))
              str = 'number'
              ierr=pio_put_att (tape(t)%File, cosp_scol_var, 'units', trim(str))

           case('cosp_htmisr')
              ierr=pio_def_var (tape(t)%File, 'cosp_htmisr', PIO_DOUBLE,  mdimids(i:i),cosp_htmisr_var)
              str = 'COSP MISR height'
              ierr=pio_put_att (tape(t)%File, cosp_htmisr_var, 'long_name', trim(str))
              str = 'km'
              ierr=pio_put_att (tape(t)%File, cosp_htmisr_var, 'units', trim(str))
              str = 'cosp_htmisr_bnds'
              ierr=pio_put_att (tape(t)%File, cosp_htmisr_var, 'bounds', trim(str))
              ierr=pio_def_var (tape(t)%File, str, PIO_DOUBLE, (/bnddim,mdimids(i:i)/), cosp_htmisr_bnds)
           case('cosp_sza')
              ierr=pio_def_var (tape(t)%File, 'cosp_sza', PIO_DOUBLE,  mdimids(i:i),cosp_sza_var)
              str = 'COSP Parasol SZA'
              ierr=pio_put_att (tape(t)%File, cosp_sza_var, 'long_name', trim(str))
              str = 'degrees'
              ierr=pio_put_att (tape(t)%File, cosp_sza_var, 'units', trim(str))
              !           case('cosp_nbnds')
              !          case('cosp_nhtml')
           case('lev','ilev')   ! just prevent the message for these cases.
           case default
              if(masterproc) write(iulog,*) 'Warning: no dimvar found for registered mdim named: ',hist_mdims(i)%name
           end select
        end if
     end do

     ! There should be no further calls to pio_def_dim after the registeredmdims loop

     ierr=pio_def_var (tape(t)%File,'P0',pio_double,ps0var)
     str = 'reference pressure'
     ierr=pio_put_att (tape(t)%File, ps0var, 'long_name', trim(str))
     ierr=pio_put_att (tape(t)%File, ps0var, 'units', 'Pa')

     call get_ref_date(yr, mon, day, nbsec)
     nbdate = yr*10000 + mon*100 + day
     ierr=pio_def_var (tape(t)%File,'time',pio_double,(/TIMDIM/),tape(t)%timeid)
     ierr=pio_put_att (tape(t)%File, tape(t)%timeid, 'long_name', 'time')
     str = 'days since ' // date2yyyymmdd(nbdate) // ' ' // sec2hms(nbsec)
     ierr=pio_put_att (tape(t)%File, tape(t)%timeid, 'units', trim(str))

     calendar = timemgr_get_calendar_cf()
     ierr=pio_put_att (tape(t)%File, tape(t)%timeid, 'calendar', calendar)

     
     ierr=pio_def_var (tape(t)%File,'date    ',pio_int,(/timdim/),tape(t)%dateid)
     str = 'current date (YYYYMMDD)'
     ierr=pio_put_att (tape(t)%File, tape(t)%dateid, 'long_name', trim(str))


     ierr=pio_def_var (tape(t)%File,'datesec ',pio_int,(/timdim/), tape(t)%datesecid)
     str = 'current seconds of current date'
     ierr=pio_put_att (tape(t)%File, tape(t)%datesecid, 'long_name', trim(str))


     !     
     ! Character header information 
     !
     str = 'CF-1.0'
     ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'Conventions', trim(str))
     ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'source', 'CAM')
#if ( defined BFB_CAM_SCAM_IOP )
     ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'CAM_GENERATED_FORCING','create SCAM IOP dataset')
#endif
     ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'case',caseid)
     ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'title',ctitle)
     ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'logname',logname)
     ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'host', host)
     ierr= pio_put_att (tape(t)%File, PIO_GLOBAL, 'Version', &
          '$Name$')
     ierr= pio_put_att (tape(t)%File, PIO_GLOBAL, 'revision_Id', &
          '$Id$')
     ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'initial_file', ncdata)
     ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'topography_file', bnd_topo)

     if(.not. is_satfile(t)) then
        if(ngroup(t)==0 .or. restart) then
           if(ncoloutput) then
              ierr=pio_def_var (tape(t)%File,'lat',pio_double,(/NCOLDIM/),latvar)
              if(dycore_is('UNSTRUCTURED')) then
                 ierr=pio_def_var (tape(t)%File,'area',pio_double,(/NCOLDIM/),areavar)
              end if
           else      
              ierr=pio_def_var (tape(t)%File,'lat',pio_double,(/LATDIM/),latvar)
           end if
           ierr=pio_put_att (tape(t)%File, latvar, 'long_name', 'latitude')
           ierr=pio_put_att (tape(t)%File, latvar, 'units', 'degrees_north')
           
           if(ncoloutput ) then
              ierr=pio_def_var (tape(t)%File,'lon',pio_double,(/NCOLDIM/),lonvar)
           else
              ierr=pio_def_var (tape(t)%File,'lon',pio_double,(/LONDIM/),lonvar)
           end if
           ierr=pio_put_att (tape(t)%File, lonvar,'long_name','longitude')
           ierr=pio_put_att (tape(t)%File, lonvar,'units','degrees_east')

           ! If a staggered grid is in use, output the lat and lon arrays variables.
           if(dycore_is('LR') .and. .not. ncoloutput) then
              
              ierr=pio_def_var (tape(t)%File, 'slat', pio_double,(/SLATDIM/),slatvar)
              ierr=pio_put_att (tape(t)%File, slatvar, 'long_name', 'staggered latitude')
              ierr=pio_put_att (tape(t)%File, slatvar, 'units', 'degrees_north')

              ierr=pio_def_var (tape(t)%File,'slon',pio_double,(/SLONDIM/),slonvar)
              ierr=pio_put_att (tape(t)%File, slonvar, 'long_name', 'staggered longitude')
              ierr=pio_put_att(tape(t)%File, slonvar, 'units', 'degrees_east')

              ierr=pio_def_var (tape(t)%File,'w_stag',pio_double,(/slatdim/),wsid)
              str = 'staggered latitude weights'
              ierr=pio_put_att (tape(t)%File, wsid, 'long_name', trim(str))
           endif
        end if

        ierr=pio_put_att (tape(t)%File, tape(t)%timeid, 'bounds', 'time_bnds')

        ierr=pio_def_var (tape(t)%File,'time_bnds',pio_double,(/bnddim,timdim/),tape(t)%tbndid)
        ierr=pio_put_att (tape(t)%File, tape(t)%tbndid, 'long_name', 'time interval endpoints')
        !
        ! Character
        !
        dimenchar(1) = chardim
        dimenchar(2) = timdim
        ierr=pio_def_var (tape(t)%File,'date_written',PIO_CHAR,dimenchar, tape(t)%date_writtenid)
        ierr=pio_def_var (tape(t)%File,'time_written',PIO_CHAR,dimenchar, tape(t)%time_writtenid)
        !
        ! Integer Header
        !
        ierr=pio_def_var (tape(t)%File,'ntrm',PIO_INT,ntrmid)
        str = 'spectral truncation parameter M'
        ierr=pio_put_att (tape(t)%File, ntrmid, 'long_name', trim(str))

        ierr=pio_def_var (tape(t)%File,'ntrn',PIO_INT, ntrnid)
        str = 'spectral truncation parameter N'
        ierr=pio_put_att (tape(t)%File, ntrnid, 'long_name', trim(str))

        ierr=pio_def_var (tape(t)%File,'ntrk',PIO_INT,ntrkid)
        str = 'spectral truncation parameter K'
        ierr=pio_put_att (tape(t)%File, ntrkid, 'long_name', trim(str))

        ierr=pio_def_var (tape(t)%File,'ndbase',PIO_INT,tape(t)%ndbaseid)
        str = 'base day'
        ierr=pio_put_att (tape(t)%File, tape(t)%ndbaseid, 'long_name', trim(str))

        ierr=pio_def_var (tape(t)%File,'nsbase',PIO_INT,tape(t)%nsbaseid)
        str = 'seconds of base day'
        ierr=pio_put_att (tape(t)%File, tape(t)%nsbaseid, 'long_name', trim(str))

        ierr=pio_def_var (tape(t)%File,'nbdate',PIO_INT,tape(t)%nbdateid)
        str = 'base date (YYYYMMDD)'
        ierr=pio_put_att (tape(t)%File, tape(t)%nbdateid, 'long_name', trim(str))

#if ( defined BFB_CAM_SCAM_IOP )
        ierr=pio_def_var (tape(t)%File,'bdate',PIO_INT,tape(t)%bdateid)
        str = 'base date (YYYYMMDD)'
        ierr=pio_put_att (tape(t)%File, tape(t)%bdateid, 'long_name', trim(str))
#endif
        ierr=pio_def_var (tape(t)%File,'nbsec',PIO_INT,tape(t)%nbsecid)
        str = 'seconds of base date'
        ierr=pio_put_att (tape(t)%File, tape(t)%nbsecid, 'long_name', trim(str))

        ierr=pio_def_var (tape(t)%File,'mdt',PIO_INT,tape(t)%mdtid)
        ierr=pio_put_att (tape(t)%File, tape(t)%mdtid, 'long_name', 'timestep')
        ierr=pio_put_att (tape(t)%File, tape(t)%mdtid, 'units', 's')


        if(latdim > -1) then
           if(.not. is_initfile(file_index=t) .and. .not. dycore_is('UNSTRUCTURED') ) then
              ierr=pio_def_var (tape(t)%File,'nlon',PIO_INT,(/latdim/),tape(t)%nlonid)
              str = 'number of longitudes'
              ierr=pio_put_att (tape(t)%File, tape(t)%nlonid, 'long_name', trim(str))

              ierr=pio_def_var (tape(t)%File,'wnummax',PIO_INT,(/latdim/),tape(t)%wnummaxid)
              str = 'cutoff Fourier wavenumber'
              ierr=pio_put_att (tape(t)%File, tape(t)%wnummaxid, 'long_name', trim(str))
           end if
        !
        ! Floating point time-invariant
        !
           ierr=pio_def_var (tape(t)%File,'gw',PIO_DOUBLE,(/latdim/),gwid)
           str = 'gauss weights'
           ierr=pio_put_att (tape(t)%File, gwid, 'long_name', trim(str))
        end if

        !
        ! Create variables for model timing and header information 
        !

        ierr=pio_def_var (tape(t)%File,'ndcur   ',pio_int,(/timdim/),tape(t)%ndcurid)
        str = 'current day (from base day)'
        ierr=pio_put_att (tape(t)%File, tape(t)%ndcurid, 'long_name', trim(str))

        ierr=pio_def_var (tape(t)%File,'nscur   ',pio_int,(/timdim/),tape(t)%nscurid)
        str = 'current seconds of current day'
        ierr=pio_put_att (tape(t)%File, tape(t)%nscurid, 'long_name', trim(str))


        if (.not. is_initfile(file_index=t)) then
           ! Don't write the GHG/Solar forcing data to the IC file.
           ierr=pio_def_var (tape(t)%File,'co2vmr  ',pio_double,(/timdim/),tape(t)%co2vmrid)
           str = 'co2 volume mixing ratio'
           ierr=pio_put_att (tape(t)%File, tape(t)%co2vmrid, 'long_name', trim(str))

           ierr=pio_def_var (tape(t)%File,'ch4vmr  ',pio_double,(/timdim/),tape(t)%ch4vmrid)
           str = 'ch4 volume mixing ratio'
           ierr=pio_put_att (tape(t)%File, tape(t)%ch4vmrid, 'long_name', trim(str))

           ierr=pio_def_var (tape(t)%File,'n2ovmr  ',pio_double,(/timdim/),tape(t)%n2ovmrid)
           str = 'n2o volume mixing ratio'
           ierr=pio_put_att (tape(t)%File, tape(t)%n2ovmrid, 'long_name', trim(str))

           ierr=pio_def_var (tape(t)%File,'f11vmr  ',pio_double,(/timdim/),tape(t)%f11vmrid)
           str = 'f11 volume mixing ratio'
           ierr=pio_put_att (tape(t)%File, tape(t)%f11vmrid, 'long_name', trim(str))

           ierr=pio_def_var (tape(t)%File,'f12vmr  ',pio_double,(/timdim/),tape(t)%f12vmrid)
           str = 'f12 volume mixing ratio'
           ierr=pio_put_att (tape(t)%File, tape(t)%f12vmrid, 'long_name', trim(str))

           ierr=pio_def_var (tape(t)%File,'sol_tsi ',pio_double,(/timdim/),tape(t)%sol_tsiid)
           str = 'total solar irradiance'
           ierr=pio_put_att (tape(t)%File, tape(t)%sol_tsiid, 'long_name', trim(str))
           str = 'W/m2'
           ierr=pio_put_att (tape(t)%File, tape(t)%sol_tsiid, 'units', trim(str))
        end if


#if ( defined BFB_CAM_SCAM_IOP )
        ierr=pio_def_var (tape(t)%File,'tsec ',pio_int,(/timdim/), tape(t)%tsecid)
        str = 'current seconds of current date needed for scam'
        ierr=pio_put_att (tape(t)%File, tape(t)%tsecid, 'long_name', trim(str))
#endif
        ierr=pio_def_var (tape(t)%File,'nsteph  ',pio_int,(/timdim/),tape(t)%nstephid)
        str = 'current timestep'
        ierr=pio_put_att (tape(t)%File, tape(t)%nstephid, 'long_name', trim(str))
     end if



     !
     ! Create variables and attributes for field list
     !

     do f=1,nflds(t)
        select case (tape(t)%hlist(f)%avgflag)
        case ('A')
           tape(t)%hlist(f)%time_op(:) = 'mean'
        case ('B')
           tape(t)%hlist(f)%time_op(:) = 'mean00z'
        case ('I')
           tape(t)%hlist(f)%time_op(:) = ' '
        case ('X')
           tape(t)%hlist(f)%time_op(:) = 'maximum'
        case ('M')
           tape(t)%hlist(f)%time_op(:) = 'minimum'
        case('L')
           tape(t)%hlist(f)%time_op(:) = LT_DESC
        case default
           call endrun ('h_define: unknown avgflag='//tape(t)%hlist(f)%avgflag)
        end select

        numlev = tape(t)%hlist(f)%field%numlev

        if (tape(t)%hlist(f)%hwrt_prec == 8) then
           ncreal = pio_double
        else
           ncreal = pio_real 
        end if
        if(.not.associated(tape(t)%hlist(f)%varid)) then
           allocate(tape(t)%hlist(f)%varid(1:max(1,ngroup(t))))
        end if
        fname_tmp = strip_suffix(tape(t)%hlist(f)%field%name)

        if(is_satfile(t)) then
           hdims=0
           nfils(t)=1
           call sat_hist_define(tape(t)%File)
        else
           if(ncoloutput) then
              dimindex(1)=ncoldim
              hdims=1
              nacsdims(1) = ncoldim
           else
              nacsdims=(/londim,latdim/)
              dimindex(1)=londim               
              dimindex(2)=latdim               
              hdims=2
              if(dycore_is('LR') .and. numlev>1) then
                 if(fname_tmp .eq. 'US' .or. fname_tmp .eq. 'FU_S') then
                    dimindex(1:2) = (/londim,slatdim/)
                    nacsdims=(/londim,slatdim/)
                 else if(fname_tmp .eq. 'VS' .or. fname_tmp .eq. 'FV_S') then
                    dimindex(1:2) = (/slondim,latdim/)
                    nacsdims=(/slondim,latdim/)
                 else
                    dimindex(1:2) = (/londim,latdim/)
                    nacsdims=(/londim,latdim/)
                 end if
              end if
           end if
        endif

        !
        !  Create variables and atributes for fields written out as columns
        !         
        if(associated(tape(t)%hlist(f)%field%mdims)) then
           mdims => tape(t)%hlist(f)%field%mdims
           mdimsize = size(mdims)
        else if(tape(t)%hlist(f)%field%numlev>1) then
           call endrun('mdims not defined for variable '//trim(tape(t)%hlist(f)%field%name))
        else
           mdimsize=0   
        endif
        do i=0,ngroup(t)
           if(i==0 .and. ngroup(t)>0 .and. .not. restart) cycle
           fname_tmp = strip_suffix(tape(t)%hlist(f)%field%name)
           if(i==0) then  ! i==0 full field output
              if(.not.associated(tape(t)%hlist(f)%varid)) then
                 allocate(tape(t)%hlist(f)%varid(1))
              end if
              varid => tape(t)%hlist(f)%varid(1)
              if(restart) then
                 if (tape(t)%hlist(f)%field%flag_xyfill .or. (avgflag_pertape(t)=='L')) then
                    ierr=pio_def_var(tape(t)%File,trim(fname_tmp)//'_nacs',pio_int,nacsdims(1:hdims),&
                         tape(t)%hlist(f)%nacs_varid)
                 else
                    ierr=pio_def_var(tape(t)%File,trim(fname_tmp)//'_nacs',pio_int,&
                         tape(t)%hlist(f)%nacs_varid)
                 end if
              endif
           else   ! i>0 column group output
              if(restart) exit
              ff=(f-1)*ngroup(t)+i
              dimindex(1)=grouplondim(i)               
              dimindex(2)=grouplatdim(i)               
              hdims=2
              if(dycore_is('LR')) then
                 if(fname_tmp .eq. 'US' .or. fname_tmp .eq. 'FU_S') then
                    dimindex(2)=grouplatdim_st(i)
                 else if(fname_tmp .eq. 'VS' .or. fname_tmp .eq. 'FV_S') then
                    dimindex(1)=grouplondim_st(i)               
                 end if
              endif
              fname_tmp = tape(t)%hlist(f)%field_column_name(i)
              varid => tape(t)%hlist(f)%varid(i)
           endif
           if(mdimsize==0) then
              if(restart) then
                 tdims=hdims
              else
                 dimindex(hdims+1) = timdim
                 tdims=hdims+1
              end if
              ierr=pio_def_var(tape(t)%File,trim(fname_tmp),ncreal,dimindex(1:tdims),varid)
           else 
              do j=1,mdimsize
                 dimindex(hdims+j) = mdimids(mdims(j))
              end do
              if(restart) then
                 tdims=hdims+mdimsize
              else
                 tdims=hdims+mdimsize+1
                 dimindex(tdims)=timdim
              end if
              ierr=pio_def_var(tape(t)%File,trim(fname_tmp),ncreal,dimindex(1:tdims),varid)
              ierr=pio_put_att(tape(t)%FIle,varid,'mdims', mdims(1:mdimsize))
           end if

           str = tape(t)%hlist(f)%field%sampling_seq
           if (len_trim(str)>0) then
              ierr=pio_put_att (tape(t)%File, varid,'Sampling_Sequence', trim(str))
           end if


           if (tape(t)%hlist(f)%field%flag_xyfill) then
              ! Add both _FillValue and missing_value to cover expectations of various applications.
              ! The attribute type must match the data type.
              if (tape(t)%hlist(f)%hwrt_prec == 8) then
                 ierr=pio_put_att (tape(t)%File, varid, '_FillValue', tape(t)%hlist(f)%field%fillvalue)
                 ierr=pio_put_att (tape(t)%File, varid, 'missing_value', tape(t)%hlist(f)%field%fillvalue)
              else
                 ierr=pio_put_att (tape(t)%File, varid, '_FillValue', real(tape(t)%hlist(f)%field%fillvalue,r4))
                 ierr=pio_put_att (tape(t)%File, varid, 'missing_value', real(tape(t)%hlist(f)%field%fillvalue,r4))
              end if
           end if

           str = tape(t)%hlist(f)%field%units
           ierr=pio_put_att (tape(t)%File, varid, 'units', trim(str))

           str = tape(t)%hlist(f)%field%long_name
           ierr=pio_put_att (tape(t)%File, varid, 'long_name', trim(str))
           !
           ! Assign field attributes defining valid levels and averaging info
           !
           str = tape(t)%hlist(f)%time_op
           select case (str)
           case ('mean', 'maximum', 'minimum' )
              ierr=pio_put_att (tape(t)%File, varid,'cell_methods', 'time: '//str)
           end select
           if(i>0) ierr= pio_put_att(tape(t)%File, varid, 'basename',tape(t)%hlist(f)%field%name)
        end do
     end do
     !
     ret = pio_enddef(tape(t)%File)

     if(masterproc) then
        write(iulog,*)'H_DEFINE: Successfully opened netcdf file '
     endif
     !
     ! Write time-invariant portion of history header
     !
     ierr = pio_put_var (tape(t)%File, ps0var, (/ps0/))
     if(.not. is_satfile(t)) then
        londeg => get_dyn_grid_parm_real2d('londeg')
        latdeg => get_dyn_grid_parm_real1d('latdeg')
        if(dycore_is('LR')) then
           londeg_st => get_dyn_grid_parm_real2d('londeg_st')
           latdeg_st => get_dyn_grid_parm_real1d('latdeg_st')
        end if
        if(ncoloutput) then
           if(dycore_is('UNSTRUCTURED')) then
              ! save memory for these global arrays by doing them 1 by 1:
              allocate(harea(ncol))
              call get_horiz_grid_d(ncol, clon_d_out=harea)
              do j=1,ncol
                 harea(j) = harea(j)*radtodeg
              end do
              ierr = pio_put_var (tape(t)%File, lonvar, harea)

              call get_horiz_grid_d(ncol, area_d_out=harea)
              ierr = pio_put_var (tape(t)%File, areavar, harea)

              call get_horiz_grid_d(ncol, clat_d_out=harea)
              do j=1,ncol
                 harea(j) = harea(j)*radtodeg
              end do
              ierr = pio_put_var (tape(t)%File, latvar, harea)
              deallocate(harea)
           else
              allocate(latdeg2d(pcols, begchunk:endchunk),londeg(pcols,begchunk:endchunk))
              latdeg2d=0_r8
              londeg=0_r8
              do lchnk=begchunk,endchunk
                 call get_rlat_all_p(lchnk,pcols,latdeg2d(:,lchnk))
                 call get_rlon_all_p(lchnk,pcols,londeg(:,lchnk))
              end do
              latdeg2d = latdeg2d*radtodeg
              londeg = londeg*radtodeg
              call get_phys_decomp(iodesc, 1,(/1/),1,pio_double,column=tape(t)%column(1))
              call pio_write_darray(tape(t)%File, lonvar, iodesc, londeg, ierr)
              call pio_write_darray(tape(t)%File, latvar, iodesc, latdeg2d, ierr)
              deallocate(latdeg2d,londeg)
           end if
        else if(interpolate) then
           ierr = pio_put_var (tape(t)%File, latvar, get_interp_lat())
           ierr = pio_put_var (tape(t)%File, lonvar, get_interp_lon())
        else if(ngroup(t)==0 .or. restart) then

           ierr=pio_put_var (tape(t)%File, latvar, latdeg)
           plon = get_dyn_grid_parm('plon')
           allocate(alon(plon))
           do i=1,plon
              alon(i) = (i-1) * 360.0_r8 / plon
           end do

           if (single_column) alon(1)=scmlon
           ierr = pio_put_var (tape(t)%File, lonvar, alon)
           deallocate(alon)
           if(dycore_is('LR')) then
              latdeg_st => get_dyn_grid_parm_real1d('latdeg_st')
              
              ierr = pio_put_var (tape(t)%File, slatvar, latdeg_st)
              
              londeg_st => get_dyn_grid_parm_real2d('londeg_st')
              ierr = pio_put_var (tape(t)%File, slonvar, londeg_st(:,1))
              
              w_staggered => get_dyn_grid_parm_real1d('w_staggered')
              ierr = pio_put_var (tape(t)%File, wsid, w_staggered)
           endif
        end if
        !
        ! write out coordinate variables for columns
        !
        if (.not. restart .and. ngroup(t)>0) then
           do i = 1, ngroup(t)
              ierr = pio_put_var (tape(t)%File, glonvar(i), &
                   londeg(tape(t)%column(i)%columnlon(1):tape(t)%column(i)%columnlon(2),1))
              ierr = pio_put_var (tape(t)%File, glatvar(i), &
                   latdeg(tape(t)%column(i)%columnlat(1):tape(t)%column(i)%columnlat(2)))
           end do
           deallocate(glonvar, glatvar,grouplatdim,grouplondim)
           !
           ! For staggered (FV) fields, write out staggered coordinate variables for columns
           !
           do f=1,nflds(t)
              select case (tape(t)%hlist(f)%field%name)
              case ('US','FU_S')
                 do i = 1, ngroup(t)
                    ierr = pio_put_var(tape(t)%File, glatvar_st(i), &
                         latdeg_st(tape(t)%column_st(i)%columnlat(1):tape(t)%column_st(i)%columnlat(2)))
                 end do
              case ('VS','FV_S')
                 do i = 1, ngroup(t)
                    ierr = pio_put_var(tape(t)%File, glonvar_st(i), &
                         londeg_st(tape(t)%column_st(i)%columnlon(1):tape(t)%column_st(i)%columnlon(2),1))
                 end do
              end select
           end do
        end if


        if(latdim>-1) then
           w => get_dyn_grid_parm_real1d('w')
           ierr = pio_put_var (tape(t)%File, gwid, w)
        end if
        ierr = pio_put_var (tape(t)%File, ntrmid, (/ptrm/))
        ierr = pio_put_var (tape(t)%File, ntrnid, (/ptrn/))
        ierr = pio_put_var (tape(t)%File, ntrkid, (/ptrk/))
        dtime = get_step_size()
        ierr = pio_put_var (tape(t)%File, tape(t)%mdtid, (/dtime/))
        !
        ! Model date info
        !
        ierr = pio_put_var (tape(t)%File, tape(t)%ndbaseid, (/ndbase/))
        ierr = pio_put_var (tape(t)%File, tape(t)%nsbaseid, (/nsbase/))
        
        ierr = pio_put_var (tape(t)%File, tape(t)%nbdateid, (/nbdate/))
#if ( defined BFB_CAM_SCAM_IOP )
        ierr = pio_put_var (tape(t)%File, tape(t)%bdateid, (/nbdate/))
#endif
        ierr = pio_put_var (tape(t)%File, tape(t)%nbsecid, (/nbsec/))
        !
        ! Reduced grid info
        !

        if(.not. is_initfile(file_index=t) .and. .not. dycore_is('UNSTRUCTURED') .and. latdim>-1) then
           ierr = pio_put_var (tape(t)%File, tape(t)%nlonid, nlon)
           ierr = pio_put_var (tape(t)%File, tape(t)%wnummaxid, wnummax)
        end if
     end if
     if(allocated(grouplondim_st)) then
        deallocate(glonvar_st, grouplondim_st)
     end if
     if(allocated(grouplatdim_st)) then
        deallocate(glatvar_st, grouplatdim_st)
     end if
     !
     ! 0.01 converts Pascals to millibars
     !
     alev(:pver) = 0.01_r8*ps0*(hyam(:pver) + hybm(:pver))
     ailev(:pverp) = 0.01_r8*ps0*(hyai(:pverp) + hybi(:pverp))


     do i=1,registeredmdims
        select case(hist_mdims(i)%name)
        case ('lev')
           ierr = pio_put_var (tape(t)%File, levvar, alev)
           ierr = pio_put_var (tape(t)%File, hyamid, hyam)
           ierr = pio_put_var (tape(t)%File, hybmid, hybm)
        case ('ilev')
           ierr = pio_put_var (tape(t)%File, ilevvar, ailev)
           ierr = pio_put_var (tape(t)%File, hyaiid, hyai)
           ierr = pio_put_var (tape(t)%File, hybiid, hybi)
        end select
        if(.not. (is_satfile(t) .or. is_initfile(t))) then
           select case(hist_mdims(i)%name)
           case('cosp_prs')
              ierr = pio_put_var (tape(t)%file, cosp_prs_var, prsmid_cosp)
              ierr = pio_put_var (tape(t)%file, cosp_prs_bnds, prslim_cosp)
           case('cosp_tau')
              ierr = pio_put_var (tape(t)%file, cosp_tau_var, taumid_cosp)
              ierr = pio_put_var (tape(t)%file, cosp_tau_bnds, taulim_cosp)
           case('cosp_tau_modis')
              ierr = pio_put_var (tape(t)%file, cosp_tau_modis_var, taumid_cosp_modis)
              ierr = pio_put_var (tape(t)%file, cosp_tau_modis_bnds, taulim_cosp_modis)
           case('cosp_ht')
              ierr = pio_put_var (tape(t)%file, cosp_ht_var, htmid_cosp)
              ierr = pio_put_var (tape(t)%file, cosp_ht_bnds, htlim_cosp)
           case('cosp_dbze')
              ierr = pio_put_var (tape(t)%file, cosp_dbze_var, dbzemid_cosp)
              ierr = pio_put_var (tape(t)%file, cosp_dbze_bnds, dbzelim_cosp)
           case('cosp_sr')
              ierr = pio_put_var (tape(t)%file, cosp_sr_var, srmid_cosp)
              ierr = pio_put_var (tape(t)%file, cosp_sr_bnds, srlim_cosp)
           case('cosp_scol')
              ierr = pio_put_var (tape(t)%file, cosp_scol_var, scol_cosp)
           case('cosp_htmisr')              
              ierr = pio_put_var (tape(t)%file, cosp_htmisr_var, htmisrmid_cosp)
              ierr = pio_put_var (tape(t)%file, cosp_htmisr_bnds, htmisrlim_cosp)
           case('cosp_sza')
              ierr = pio_put_var (tape(t)%file, cosp_sza_var, sza_cosp)
           case('lev','ilev')  ! avoids the warning message
           case default
              if(masterproc) write(iulog,*) 'Warning: no dimvar found for registered mdim named: ',hist_mdims(i)%name
           end select
        end if
     end do




   end subroutine h_define

!#######################################################################


   subroutine h_normalize (f, t)

   use dycore, only: dycore_is

!
!----------------------------------------------------------------------- 
! 
! Purpose: Normalize fields on a history file by the number of accumulations
! 
! Method: Loop over fields on the tape.  Need averaging flag and number of
!         accumulations to perform normalization.
! 
!-----------------------------------------------------------------------
!
! Input arguments
!
      integer, intent(in) :: f       ! field index
      integer, intent(in) :: t       ! tape index
!
! Local workspace
!
      integer begdim2                ! on-node vert start index
      integer enddim2                ! on-node vert end index
      integer c                      ! chunk (or lat) index
      integer endi                   ! terminating column index
      integer begdim1                ! on-node dim1 start index
      integer enddim1                ! on-node dim1 end index  
      integer begdim3                ! on-node chunk or lat start index
      integer enddim3                ! on-node chunk or lat end index  
      integer :: ib, ie    ! beginning and ending indices of first dimension
      integer :: jb, je    ! beginning and ending indices of second dimension
      integer k

      logical :: flag_xyfill         ! non-applicable xy points flagged with fillvalue
      character*1 avgflag            ! averaging flag

      call t_startf ('h_normalize')

      begdim2 = tape(t)%hlist(f)%field%begdim2
      enddim2 = tape(t)%hlist(f)%field%enddim2
      avgflag = tape(t)%hlist(f)%avgflag
!
! normalize by number of accumulations for averaged case
!
      begdim1 = tape(t)%hlist(f)%field%begdim1
      enddim1 = tape(t)%hlist(f)%field%enddim1
      begdim3 = tape(t)%hlist(f)%field%begdim3
      enddim3 = tape(t)%hlist(f)%field%enddim3
      flag_xyfill = tape(t)%hlist(f)%field%flag_xyfill
      do c=begdim3,enddim3
 ! Some slices might not be entirely full  -- use endi instead of enddim2
         if(associated(tape(t)%hlist(f)%field%colperdim3)) then         
            endi    = tape(t)%hlist(f)%field%colperdim3(c)
         else
            endi = enddim1-begdim1+1 
         end if

         ib = begdim1
         ie = begdim1+endi-1
         jb = begdim2
         je = enddim2
         
         if (flag_xyfill) then
               do k=jb,je
                  where (tape(t)%hlist(f)%nacs(ib:ie,c) == 0)
                     tape(t)%hlist(f)%hbuf(ib:ie,k,c) = tape(t)%hlist(f)%field%fillvalue
                  endwhere
               end do
         end if

         if (avgflag == 'A' .or. avgflag == 'B') then
            if(flag_xyfill) then
               do k=jb,je
                  where (tape(t)%hlist(f)%nacs(ib:ie,c) /= 0)
                     tape(t)%hlist(f)%hbuf(ib:ie,k,c) = &
                          tape(t)%hlist(f)%hbuf(ib:ie,k,c) &
                          / tape(t)%hlist(f)%nacs(ib:ie,c)
                  endwhere
               end do
            else if(tape(t)%hlist(f)%nacs(1,c) > 0) then
               do k=jb,je
                  tape(t)%hlist(f)%hbuf(ib:ie,k,c) = &
                       tape(t)%hlist(f)%hbuf(ib:ie,k,c) &
                       / tape(t)%hlist(f)%nacs(1,c)
               end do
            end if
         end if
      end do

      call t_stopf ('h_normalize')
      
      return
   end subroutine h_normalize

!#######################################################################

   subroutine h_zero (f, t)
   use dycore, only: dycore_is
!
!----------------------------------------------------------------------- 
! 
! Purpose: Zero out accumulation buffers for a tape
! 
! Method: Loop through fields on the tape
! 
!-----------------------------------------------------------------------
!
      integer, intent(in) :: f     ! field index
      integer, intent(in) :: t     ! tape index
!
! Local workspace
!
      integer begdim1              ! on-node 1st dim start index
      integer enddim1              ! on-node 1st dim end index
      integer begdim2              ! on-node vert start index
      integer enddim2              ! on-node vert end index
      integer c                    ! chunk index
      integer endi                 ! terminating column index
      integer begdim3              ! on-node chunk or lat start index
      integer enddim3              ! on-node chunk or lat end index  
      type (dim_index_3d) :: dimind ! 3-D dimension index
      
      call t_startf ('h_zero')

      begdim1 = tape(t)%hlist(f)%field%begdim1
      enddim1 = tape(t)%hlist(f)%field%enddim1
      begdim2 = tape(t)%hlist(f)%field%begdim2
      enddim2 = tape(t)%hlist(f)%field%enddim2
      begdim3 = tape(t)%hlist(f)%field%begdim3
      enddim3 = tape(t)%hlist(f)%field%enddim3

      if(associated(tape(t)%hlist(f)%field%colperdim3)) then         
         do c=begdim3,enddim3
            endi    = tape(t)%hlist(f)%field%colperdim3(c)
            tape(t)%hlist(f)%hbuf(begdim1:begdim1+endi-1,begdim2:enddim2,c)=0._r8
            tape(t)%hlist(f)%nacs(:,:) = 0
         end do
      else
         endi = tape(t)%hlist(f)%field%enddim1 - tape(t)%hlist(f)%field%begdim1 + 1
         tape(t)%hlist(f)%hbuf(begdim1:begdim1+endi-1,begdim2:enddim2,begdim3:enddim3)=0._r8
         tape(t)%hlist(f)%nacs(:,:) = 0
      end if

      call t_stopf ('h_zero')

      return
   end subroutine h_zero

!#######################################################################

   subroutine dump_field (f, t, restart)
     use interp_mod,    only: latlon_interpolation
     use cam_pio_utils, only: get_decomp, pio_subsystem     
     use dycore,        only: dycore_is
     use interp_mod, only : write_interpolated
     use spmd_utils, only : iam
     integer, intent(in) :: f, t
     logical, intent(in) :: restart
!
!----------------------------------------------------------------------- 
! 
! Purpose: Write a variable to a history tape
! 
! Method: If SPMD, first gather the data to the master processor.  
!         Next, transpose the data to COORDS order (the default).
!         Finally, issue the netcdf call to write the variable
! 
!-----------------------------------------------------------------------
     integer :: bsize, ierr
     
     real(r4), allocatable :: localbuf4(:)
     type(io_desc_t), pointer :: int_iodesc, iodesc
     type(var_desc_t), pointer :: varid, varid_vecv
     character(len=max_fieldname_len) :: fname
     type(column_info) :: tmpcolumn
     integer :: ff

     interpolate = (latlon_interpolation(t) .and. .not. restart )
     do i=1,max(ngroup(t),1)
        if(restart .and. i>1) exit
        varid=>tape(t)%hlist(f)%varid(i)

        if(restart) then
           call pio_setframe(varid, int(1,kind=PIO_Offset))
           if(.not. interpolate) then
              call get_decomp(iodesc, tape(t)%hlist(f)%field, pio_double)
           end if
        else
           call pio_setframe(varid, int(max(1,nfils(t)),kind=PIO_Offset))
           if(ngroup(t)>0) then
              tmpcolumn = tape(t)%column(i)
              nullify(tmpcolumn%hmask)
              if(dycore_is('LR')) then
                 fname = tape(t)%hlist(f)%field%name
                 if(fname .eq. 'US' .or. fname .eq. 'FU_S') then
                    tmpcolumn%columnlat(1) = tape(t)%column_st(i)%columnlat(1)+1 ! Adding 1 because staggered lat indexing for "xzybuf" starts at 2
                    tmpcolumn%columnlat(2) = tape(t)%column_st(i)%columnlat(2)+1
                    tmpcolumn%num_lats =  tape(t)%column_st(i)%num_lats
                 else if(fname .eq. 'VS' .or. fname .eq. 'FV_S') then
                    tmpcolumn%columnlon = tape(t)%column_st(i)%columnlon
                    tmpcolumn%num_lons =  tape(t)%column_st(i)%num_lons
                 end if
              end if
              allocate(iodesc)  ! no lookup table for columns
              call get_decomp(iodesc, tape(t)%hlist(f)%field, pio_double,column=tmpcolumn)

#ifdef HDEBUG
              write(iulog,*)'DUMP_FIELD:',__LINE__,' writing column time indx ', &
                   nfils(t), ' field id', tape(t)%hlist(f)%varid(i)%varid, &
                   ' name ', trim(tape(t)%hlist(f)%field_column_name(i)), &
                   'numlons=',tmpcolumn%num_lons,'numlats=', tmpcolumn%num_lats, ' base name ', &
                   trim(tape(t)%hlist(f)%field%name),' columnlon=', &
                   tmpcolumn%columnlon, ' columnlat=', tmpcolumn%columnlat
#endif
           else if(associated(tape(t)%column(1)%hmask)) then
              call get_decomp(iodesc, tape(t)%hlist(f)%field, pio_double, column=tape(t)%column(1))
           else if(.not. interpolate) then
              call get_decomp(iodesc, tape(t)%hlist(f)%field, pio_double)
           endif
        end if


        bsize=size(tape(t)%hlist(f)%hbuf)
        if(interpolate) then
           ff = tape(t)%hlist(f)%field%vector_compliment
           ! ff < 0 is the vector compliment of ff*(-1)
           if(ff > 0) then
              varid_vecv=>tape(t)%hlist(ff)%varid(i)
              call pio_setframe(varid_vecv, int(max(1,nfils(t)),kind=PIO_Offset))
              call write_interpolated(tape(t)%File, varid, varid_vecv, tape(t)%hlist(f)%hbuf,&
                   tape(t)%hlist(ff)%hbuf, tape(t)%hlist(f)%field%numlev,  pio_double, tape(t)%hlist(f)%field%decomp_type)
           else if(ff==0) then
              call write_interpolated(tape(t)%File, varid, tape(t)%hlist(f)%hbuf,&
                   tape(t)%hlist(f)%field%numlev,  pio_double, tape(t)%hlist(f)%field%decomp_type)
           endif
        else
           if(tape(t)%hlist(f)%field%numlev>1 .and. .not. associated(tape(t)%hlist(f)%field%mdims)) then
              call endrun('mdims not defined for variable '//trim(tape(t)%hlist(f)%field%name))
           end if
           call pio_write_darray(tape(t)%File, varid, iodesc, tape(t)%hlist(f)%hbuf, ierr)
        endif

        if(restart) then
           if(tape(t)%hlist(f)%field%flag_xyfill .or. (avgflag_pertape(t) == 'L')) then
              call get_decomp(int_iodesc, tape(t)%hlist(f)%field, pio_int, 1)
              call pio_write_darray(tape(t)%File, tape(t)%hlist(f)%nacs_varid, &
                   int_iodesc, reshape(tape(t)%hlist(f)%nacs, (/size(tape(t)%hlist(f)%nacs)/)), ierr)

           else
              ierr = pio_put_var(tape(t)%File, tape(t)%hlist(f)%nacs_varid,&
                   tape(t)%hlist(f)%nacs(1, tape(t)%hlist(f)%field%begdim3))
           end if
        end if
        if(.not. restart .and. ngroup(t)>0) then
           call pio_freedecomp(pio_subsystem, iodesc)
           deallocate(iodesc) ! no lookup table for columns
        end if

     end do
     
     return
   end subroutine dump_field

!#######################################################################

   logical function write_inithist ()
!
!-----------------------------------------------------------------------
! 
! Purpose: Set flags that will initiate dump to IC file when OUTFLD and
! WSHIST are called
! 
!-----------------------------------------------------------------------
!
      use time_manager, only: get_nstep, get_curr_date, get_step_size, is_last_step
!
! Local workspace
!
      integer :: yr, mon, day      ! year, month, and day components of
                                   ! a date
      integer :: nstep             ! current timestep number
      integer :: ncsec             ! current time of day [seconds]
      integer :: dtime             ! timestep size

!-----------------------------------------------------------------------

      write_inithist  = .false.

      if(is_initfile()) then

         nstep = get_nstep()
         call get_curr_date(yr, mon, day, ncsec)

         if    (inithist == '6-HOURLY') then
            dtime  = get_step_size()
            write_inithist = nstep /= 0 .and. mod( nstep, nint((6._r8*3600._r8)/dtime) ) == 0
         elseif(inithist == 'DAILY'   ) then
            write_inithist = nstep /= 0 .and. ncsec == 0
         elseif(inithist == 'MONTHLY' ) then
            write_inithist = nstep /= 0 .and. ncsec == 0 .and. day == 1
         elseif(inithist == 'YEARLY'  ) then
            write_inithist = nstep /= 0 .and. ncsec == 0 .and. day == 1 .and. mon == 1
         elseif(inithist == 'CAMIOP'  ) then
            write_inithist = nstep == 0 
         elseif(inithist == 'ENDOFRUN'  ) then
            write_inithist = nstep /= 0 .and. is_last_step()
         end if
      end if

      return
   end function write_inithist

!#######################################################################

   subroutine wshist (rgnht_in)
     !
     !----------------------------------------------------------------------- 
     ! 
     ! Purpose: Driver routine to write fields on history tape t
     ! 
     ! Method: For variables which do not need to be gathered (SPMD) just issue the netcdf call
     !         For those that do need to be gathered, call "dump_field" to do the operation.
     !         Finally, zero the history buffers for each field written.
     ! 
     ! Author: CCM Core Group
     ! 
     !-----------------------------------------------------------------------
     use dycore, only: dycore_is

     use time_manager, only: get_nstep, get_curr_date, get_curr_time,get_step_size
     use chem_surfvals, only: chem_surfvals_get, chem_surfvals_co2_rad
     use solar_data,    only: sol_tsi
     use sat_hist,      only: sat_hist_write
     logical, intent(in), optional :: rgnht_in(ptapes)
     !
     ! Local workspace
     !
     character(len=8) :: cdate  ! system date
     character(len=8) :: ctime  ! system time

     logical  :: rgnht(ptapes), restart
     integer t, f               ! tape, field indices
     integer ret                ! return value from netcdf call
     integer start              ! starting index required by nf_put_vara
     integer count1             ! count values required by nf_put_vara
     integer startc(2)          ! start values required by nf_put_vara (character)
     integer countc(2)          ! count values required by nf_put_vara (character)
#ifdef HDEBUG
     !      integer begdim3
     !      integer enddim3
#endif

     integer :: yr, mon, day      ! year, month, and day components of a date
     integer :: nstep             ! current timestep number
     integer :: ncdate            ! current date in integer format [yyyymmdd]
     integer :: ncsec             ! current time of day [seconds]
     integer :: ndcur             ! day component of current time
     integer :: nscur             ! seconds component of current time
     real(r8) :: time             ! current time
     real(r8) :: tdata(2)         ! time interval boundaries
     character(len=max_string_len) :: fname ! Filename
     logical :: prev              ! Label file with previous date rather than current
     integer :: ierr
     integer :: bsize
     integer :: dim1s,dim2s,ncol
#if ( defined BFB_CAM_SCAM_IOP )
     integer :: tsec             ! day component of current time
     integer :: dtime            ! seconds component of current time
#endif
     !-----------------------------------------------------------------------

     if(present(rgnht_in)) then
        rgnht=rgnht_in
        restart=.true.
        tape => restarthistory_tape
     else
        rgnht=.false.
        restart=.false.
        tape => history_tape
     end if

     nstep = get_nstep()
     call get_curr_date(yr, mon, day, ncsec)
     ncdate = yr*10000 + mon*100 + day
     call get_curr_time(ndcur, nscur)
     !
     ! Write time-varying portion of history file header
     !
     do t=1,ptapes
        if (nflds(t) == 0 .or. (restart .and.(.not.rgnht(t)))) cycle
        !
        ! Check if this is the IC file and if it's time to write.
        ! Else, use "nhtfrq" to determine if it's time to write
        ! the other history files.
        !
        if((.not. restart) .or. rgnht(t)) then
           if( is_initfile(file_index=t) ) then
              hstwr(t) =  write_inithist()
              prev     = .false.
           else 
              if (nhtfrq(t) == 0) then
                 hstwr(t) = nstep /= 0 .and. day == 1 .and. ncsec == 0
                 prev     = .true.
              else
                 hstwr(t) = mod(nstep,nhtfrq(t)) == 0
                 prev     = .false.
              end if
           end if
        end if
        if (hstwr(t) .or. (restart .and. rgnht(t))) then
           if(masterproc) then
              if(is_initfile(file_index=t)) then
                 write(iulog,100) yr,mon,day,ncsec
100              format('WSHIST: writing time sample to Initial Conditions h-file', &
                      ' DATE=',i4.4,'/',i2.2,'/',i2.2,' NCSEC=',i6)
              else if(is_satfile(t)) then
                 write(iulog,150) nfils(t),t,yr,mon,day,ncsec
150              format('WSHIST: writing sat columns ',i6,' to h-file ', &
                      i1,' DATE=',i4.4,'/',i2.2,'/',i2.2,' NCSEC=',i6)
              else if(hstwr(t)) then
                 write(iulog,200) nfils(t),t,yr,mon,day,ncsec
200              format('WSHIST: writing time sample ',i3,' to h-file ', &
                      i1,' DATE=',i4.4,'/',i2.2,'/',i2.2,' NCSEC=',i6)
              else if(restart .and. rgnht(t)) then
                 write(iulog,300) nfils(t),t,yr,mon,day,ncsec
300              format('WSHIST: writing history restart ',i3,' to hr-file ', &
                      i1,' DATE=',i4.4,'/',i2.2,'/',i2.2,' NCSEC=',i6)
              end if
              write(iulog,*)
           end if
           !
           ! Starting a new volume => define the metadata
           !
           if (nfils(t)==0 .or. (restart.and.rgnht(t))) then
              if(restart) then
                 fname = interpret_filename_spec( rhfilename_spec, number=(t-1))
                 hrestpath(t)=fname
              else if(is_initfile(file_index=t)) then
                 fname = interpret_filename_spec( hfilename_spec(t) )
              else
                 fname = interpret_filename_spec( hfilename_spec(t), number=(t-1), &
                      prev=prev )
              end if
              !
              ! Check that this new filename isn't the same as a previous or current filename
              !
              do f = 1, ptapes
                 if (masterproc.and. trim(fname) == trim(nhfil(f)) )then
                    write(iulog,*)'WSHIST: New filename same as old file = ', trim(fname)
                    write(iulog,*)'Is there an error in your filename specifiers?'
                    write(iulog,*)'hfilename_spec(', t, ') = ', hfilename_spec(t)
                    if ( t /= f )then
                       write(iulog,*)'hfilename_spec(', f, ') = ', hfilename_spec(f)
                    end if
                    call endrun
                 end if
              end do
              if(.not. restart) then
                 nhfil(t) = fname
                 if(masterproc) write(iulog,*)'WSHIST: nhfil(',t,')=',trim(nhfil(t))
                 cpath(t) = nhfil(t)
                 if ( len_trim(nfpath(t)) == 0 ) nfpath(t) = cpath(t)
              end if
              call h_define (t, restart)
           end if

           if(is_satfile(t)) then
              call sat_hist_write( tape(t), nflds(t), nfils(t))
           else
              if(restart) then
                 start=1
              else
                 nfils(t) = nfils(t) + 1
                 start = nfils(t)
              end if
              count1 = 1


              ierr = pio_put_var (tape(t)%File, tape(t)%ndcurid,(/start/), (/count1/),(/ndcur/))
              ierr = pio_put_var (tape(t)%File, tape(t)%nscurid,(/start/), (/count1/),(/nscur/))
              ierr = pio_put_var (tape(t)%File, tape(t)%dateid,(/start/), (/count1/),(/ncdate/))

              if (.not. is_initfile(file_index=t)) then
                 ! Don't write the GHG/Solar forcing data to the IC file.
                 ierr=pio_put_var (tape(t)%File, tape(t)%co2vmrid,(/start/), (/count1/),(/chem_surfvals_co2_rad(vmr_in=.true.)/))
                 ierr=pio_put_var (tape(t)%File, tape(t)%ch4vmrid,(/start/), (/count1/),(/chem_surfvals_get('CH4VMR')/))
                 ierr=pio_put_var (tape(t)%File, tape(t)%n2ovmrid,(/start/), (/count1/),(/chem_surfvals_get('N2OVMR')/))
                 ierr=pio_put_var (tape(t)%File, tape(t)%f11vmrid,(/start/), (/count1/),(/chem_surfvals_get('F11VMR')/))
                 ierr=pio_put_var (tape(t)%File, tape(t)%f12vmrid,(/start/), (/count1/),(/chem_surfvals_get('F12VMR')/))
                 ierr=pio_put_var (tape(t)%File, tape(t)%sol_tsiid,(/start/), (/count1/),(/sol_tsi/))
              end if



              ierr = pio_put_var (tape(t)%File, tape(t)%datesecid,(/start/),(/count1/),(/ncsec/))
#if ( defined BFB_CAM_SCAM_IOP )
              dtime = get_step_size()
              tsec=dtime*nstep
              ierr = pio_put_var (tape(t)%File, tape(t)%tsecid,(/start/),(/count1/),(/tsec/))
#endif
              ierr = pio_put_var (tape(t)%File, tape(t)%nstephid,(/start/),(/count1/),(/nstep/))
              time = ndcur + nscur/86400._r8
              ierr=pio_put_var (tape(t)%File, tape(t)%timeid, (/start/),(/count1/),(/time/))

              startc(1) = 1
              startc(2) = start
              countc(1) = 2
              countc(2) = 1
              if (is_initfile(file_index=t)) then
                 tdata = time   ! Inithist file is always instantanious data
              else
                 tdata(1) = beg_time(t)
                 tdata(2) = time
              end if
              ierr=pio_put_var (tape(t)%File, tape(t)%tbndid, startc, countc, tdata)
              if(.not.restart) beg_time(t) = time  ! update beginning time of next interval
              startc(1) = 1
              startc(2) = start
              countc(1) = 8
              countc(2) = 1
              call datetime (cdate, ctime)
              ierr = pio_put_var (tape(t)%File, tape(t)%date_writtenid, startc, countc, (/cdate/))
              ierr = pio_put_var (tape(t)%File, tape(t)%time_writtenid, startc, countc, (/ctime/))

              if(.not. restart) then
                 !$OMP PARALLEL DO PRIVATE (F)
                 do f=1,nflds(t)
                    ! Normalized averaged fields
                    if (tape(t)%hlist(f)%avgflag /= 'I') then
                       call h_normalize (f, t)
                    end if
                 end do
              end if
              !
              ! Write field to history tape.  Note that this is NOT threaded due to netcdf limitations
              !
              call t_startf ('dump_field')
              do f=1,nflds(t)
                 call dump_field(f, t, restart)
              end do
              call t_stopf ('dump_field')
              !
              ! Zero history buffers and accumulators now that the fields have been written.
              !



              if(restart) then
                 do f=1,nflds(t)
                    if(associated(tape(t)%hlist(f)%varid)) then
                       deallocate(tape(t)%hlist(f)%varid)
                       nullify(tape(t)%hlist(f)%varid)
                    end if
                 end do
                 call pio_closefile(tape(t)%File)
              else
                 !$OMP PARALLEL DO PRIVATE (F)
                 do f=1,nflds(t)
                    call h_zero (f, t)
                 end do
              end if
           end if
        end if

     end do

     return
   end subroutine wshist

!#######################################################################

   subroutine addvar (ncid, name, xtype , ndims , dimids, vid)

!
!----------------------------------------------------------------------- 
! 
! Purpose: Issue the netcdf call to add a variable to the dataset
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
!
! Input arguments
!
      type(file_desc_t), intent(inout) :: ncid       ! netcdf file id
      integer, intent(in) :: xtype      ! netcdf type flag
      integer, intent(in) :: ndims      ! number of dimensions
      integer, intent(in) :: dimids(:)  ! dimension ids
      type(var_desc_t), intent(out) :: vid        ! variable ids
      
      character(len=*), intent(in) :: name ! variable name
      integer :: ierr
!
! Local workspace
!
      type(master_entry), pointer :: listentry

      ierr=pio_def_var (ncid, name, xtype , dimids(1:ndims), vid)

      listentry=> get_entry_by_name(masterlinkedlist, name)
      if(associated(listentry)) then
         ierr=pio_put_att (ncid, vid, 'long_name', listentry%field%long_name)
         ierr=pio_put_att (ncid, vid, 'units', listentry%field%units)
      end if
         
!
! Field not found in masterlist: long_name and units attributes unknown so just
! return
!
      return
   end subroutine addvar

!#######################################################################

   subroutine addfld (fname, units, numlev, avgflag, long_name, &
                      decomp_type, flag_xyfill,sampling_seq, &
                      begdim1, enddim1, &
                      begdim2, enddim2,&
                      begdim3, enddim3, mdimnames, fill_value )
                 
!
!----------------------------------------------------------------------- 
! 
! Purpose: Add a field to the master field list
! 
! Method: Put input arguments of field name, units, number of levels, averaging flag, and 
!         long name into a type entry in the global master field list (masterlist).
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
      use ppgrid,        only: begchunk, endchunk, pver, pverp
      use rgrid,         only: nlon
      use phys_grid,     only: get_ncols_p, physgrid_set
      use dycore,        only: dycore_is
      use cam_pio_utils, only: phys_decomp, dyn_decomp, dyn_stagger_decomp
      use cam_history_support, only : fillvalue

!
! Arguments
!
      character(len=*), intent(in) :: fname      ! field name--should be "max_fieldname_len" characters long
                                                 ! or less
      character(len=*), intent(in) :: units      ! units of fname--should be 8 chars
      character(len=1), intent(in) :: avgflag    ! averaging flag
      character(len=*), intent(in) :: long_name  ! long name of field
      
      integer, intent(in) :: numlev              ! number of vertical levels (dimension and loop)
      integer, intent(in) :: decomp_type         ! decomposition type

      logical, intent(in), optional :: flag_xyfill ! non-applicable xy points flagged with fillvalue

      character(len=*), intent(in), optional :: sampling_seq ! sampling sequence - if not every timestep, 
                                                             ! how often field is sampled:  
                                                             ! every other; only during LW/SW radiation calcs, etc.
      integer, intent(in), optional :: begdim1, enddim1
      integer, intent(in), optional :: begdim2, enddim2
      integer, intent(in), optional :: begdim3, enddim3

      character(len=*), intent(in), optional :: mdimnames(:)
      real(r8), intent(in), optional :: fill_value


      integer dim1s,dim2s                        ! global size of the first and second horizontal dim.

!
! Local workspace
!
      character(len=max_fieldname_len) :: fname_tmp ! local copy of fname
      type(master_entry), pointer :: listentry

      integer :: c             ! chunk (physics ) or latitude (dynamics) index
      integer :: ncol          ! number of columns per chunk
      integer :: beglat, endlat
      integer :: beglatxy, endlatxy, beglonxy, endlonxy, beglon, endlon
      integer :: dimcnt

      beglat = get_dyn_grid_parm('beglat')
      endlat = get_dyn_grid_parm('endlat')
      beglon = get_dyn_grid_parm('beglon')
      endlon = get_dyn_grid_parm('endlon')
      beglatxy = get_dyn_grid_parm('beglatxy')
      endlatxy = get_dyn_grid_parm('endlatxy')
      beglonxy = get_dyn_grid_parm('beglonxy')
      endlonxy = get_dyn_grid_parm('endlonxy')
!
! Ensure that required grid info is available now
!
      select case (decomp_type)
      case (phys_decomp)
         if (.not. physgrid_set) then
            call endrun ('ADDFLD: Attempt to add field '//fname//' to masterlist before physics grid set')
         end if
      case (dyn_decomp)
         if (.not. dyndecomp_set) then
            call endrun ('ADDFLD: Attempt to add field '//fname//' to masterlist before dynamics grid set')
         end if
      end select
      call get_horiz_grid_dim_d(dim1s,dim2s)

!
! Ensure that new field name is not all blanks
!
      if (len_trim(fname)==0) then
         call endrun('ADDFLD: blank field name not allowed')
      end if
!
! Ensure that new field name is not longer than allowed
! (strip "&IC" suffix if it exists)
!
      fname_tmp  = fname
      fname_tmp  = strip_suffix(fname_tmp)

      if (len_trim(fname_tmp) > fieldname_len) then
         write(iulog,*)'ADDFLD: field name cannot be longer than ',fieldname_len,' characters long'
         write(iulog,*)'Field name:  ',fname
         call endrun()
      end if
!
! Ensure that new field doesn't already exist
!
      listentry => get_entry_by_name(masterlinkedlist, fname)
      if(associated(listentry)) then
         call endrun ('ADDFLD:  '//fname//' already on list')
      end if

!
! Add field to Master Field List arrays fieldn and iflds
!
      allocate(listentry)
      listentry%field%name        = fname
      listentry%field%long_name   = long_name
      listentry%field%units       = units      
      listentry%field%numlev      = numlev
      listentry%field%decomp_type = decomp_type      

      nullify(listentry%field%mdims)
      nullify(listentry%field%colperdim3)

      listentry%htapeindx(:)=-1
      listentry%act_sometape = .false.
      listentry%actflag(:) = .false.

!
! Indicate sampling sequence of field (i.e., how often "outfld" is called)
! If not every timestep (default), then give a descriptor indicating the
! sampling pattern.  Currently, the only valid value is "rad_lwsw" for sampling
! during LW/SW radiation timesteps only
!
      if (present(sampling_seq)) then
         listentry%field%sampling_seq = sampling_seq
      else
         listentry%field%sampling_seq = ' '
      end if
!
! Whether to apply xy fillvalue: default is false
!
      if (present(flag_xyfill)) then
         listentry%field%flag_xyfill = flag_xyfill
      else
         listentry%field%flag_xyfill = .false.
      end if

!
!    Allow external packages to have fillvalues different than default
!

      if(present(fill_value)) then
         listentry%field%fillvalue = fill_value
      else
         listentry%field%fillvalue = fillvalue
      endif

!
! ID and size of extra dimensions
!

      if(present(mdimnames)) then
         dimcnt = size(mdimnames)
         allocate(listentry%field%mdims(dimcnt))
         call lookup_hist_mdim_indicies(mdimnames, listentry%field%mdims)
         if(dimcnt > maxvarmdims) maxvarmdims = dimcnt
      else if(numlev>1) then
         
         allocate(listentry%field%mdims(1))
         select case(numlev)
         case(pver)
            call lookup_hist_mdim_indicies((/'lev'/),listentry%field%mdims)
         case(pverp)
            call lookup_hist_mdim_indicies((/'ilev'/),listentry%field%mdims)
         case default
            write(iulog,*) 'No Vertical dimension = ',numlev,' defined for var ',trim(fname)
            call endrun()
         end select
      endif
!
! Dimension history info based on decomposition type (dynamics or physics)
!
      nullify (listentry%field%colperdim3)
      select case (decomp_type)
      case (phys_decomp)
         listentry%field%begdim1  = 1
         listentry%field%enddim1  = pcols
         listentry%field%begdim2  = 1
         listentry%field%enddim2  = numlev
         listentry%field%begdim3  = begchunk
         listentry%field%enddim3  = endchunk
         allocate (listentry%field%colperdim3(begchunk:endchunk))
         do c=begchunk,endchunk
            ncol = get_ncols_p(c)
            listentry%field%colperdim3(c) = ncol
         end do
      case (dyn_stagger_decomp)
            listentry%field%begdim1  = beglonxy
            listentry%field%enddim1  = endlonxy
            listentry%field%begdim2  = 1
            listentry%field%enddim2  = numlev
            listentry%field%begdim3  = beglatxy
            listentry%field%enddim3  = endlatxy
      case (dyn_decomp)
         if ( dycore_is('LR') )then
            listentry%field%begdim1  = beglonxy
            listentry%field%enddim1  = endlonxy
            listentry%field%begdim2  = 1
            listentry%field%enddim2  = numlev
            listentry%field%begdim3  = beglatxy
            listentry%field%enddim3  = endlatxy
         else
            listentry%field%begdim3  = beglat
            listentry%field%enddim3  = endlat
            listentry%field%begdim2  = 1
            listentry%field%enddim2  = numlev
            listentry%field%begdim1  = 1
            listentry%field%enddim1  = dim1s
            if(.not. dycore_is('UNSTRUCTURED')) then
               allocate (listentry%field%colperdim3(beglat:endlat))
               do c=beglat,endlat
                  listentry%field%colperdim3(c) = nlon(c)
               end do
            endif
         endif

      case default
         write(iulog,*)'ADDFLD: unknown decomp_type=', decomp_type
         call endrun ()
      end select

      if(present(begdim1)) then
         listentry%field%begdim1  = begdim1
      end if
      if(present(enddim1)) then
         listentry%field%enddim1  = enddim1
      end if
      if(present(begdim2)) then
         listentry%field%begdim2  = begdim2
      end if
      if(present(enddim2)) then
         listentry%field%enddim2  = enddim2
      end if
      if(present(begdim3)) then
         listentry%field%begdim3  = begdim3
      end if
      if(present(enddim3)) then
         listentry%field%enddim3  = enddim3
      end if



!
! These 2 fields are used only in master field list, not runtime field list
!
      listentry%avgflag(:) = avgflag
      listentry%actflag(:) = .false.

      select case (avgflag)
      case ('A')
         listentry%time_op(:) = 'mean'
      case ('B')
         listentry%time_op(:) = 'mean00z'
      case ('I')
         listentry%time_op(:) = ' '
      case ('X')
         listentry%time_op(:) = 'maximum'
      case ('M')
         listentry%time_op(:) = 'minimum'
      case ('L')
         listentry%time_op(:) = LT_DESC
      case default
         call endrun ('ADDFLD: unknown avgflag='//avgflag)
      end select
      nullify(listentry%next_entry)

      call add_entry_to_master(listentry)
      return
   end subroutine addfld

   subroutine add_entry_to_master( newentry)
     type(master_entry), target, intent(in) :: newentry
     type(master_entry), pointer :: listentry
 
     if(associated(masterlinkedlist)) then
        listentry => masterlinkedlist
        do while(associated(listentry%next_entry))
           listentry=>listentry%next_entry
        end do
        listentry%next_entry=>newentry
     else
        masterlinkedlist=>newentry
     end if

   end subroutine add_entry_to_master

!#######################################################################

   subroutine wrapup (rstwr, nlend)
!
!-----------------------------------------------------------------------
!
! Purpose: 
! Close history files.
! 
! Method: 
! This routine will close any full hist. files
! or any hist. file that has data on it when restart files are being 
! written.
! If a partially full history file was disposed (for restart 
! purposes), then wrapup will open that unit back up and position 
! it for appending new data. 
!
! Original version: CCM2
!
!-----------------------------------------------------------------------
!
      use pio, only : pio_file_is_open
      use shr_kind_mod,  only: r8 => shr_kind_r8
      use pspect
      use ioFileMod
      use time_manager,  only: get_nstep, get_curr_date, get_curr_time
      use cam_pio_utils, only: cam_pio_openfile

!
! Input arguments
!
      logical, intent(in) :: rstwr   ! true => restart files are written this timestep
      logical, intent(in) :: nlend   ! Flag if time to end

!
! Local workspace
!
      integer :: nstep          ! current timestep number
      integer :: ncdate         ! current date in integer format [yyyymmdd]
      integer :: ncsec          ! time of day relative to current date [seconds]
      integer :: ndcur          ! days component of current time
      integer :: nscur          ! seconds component of current time
      integer :: yr, mon, day   ! year, month, day components of a date

      logical lfill   (ptapes)  ! Is history file ready to dispose?
      logical hdispose(ptapes)  ! Primary file disposed
      logical lhdisp            ! true => history file is disposed
      logical lhfill            ! true => history file is full

      integer t                 ! History file number
      integer f
      integer :: ierr
      real(r8) tday             ! Model day number for printout
!-----------------------------------------------------------------------

      tape => history_tape

      nstep = get_nstep()
      call get_curr_date(yr, mon, day, ncsec)
      ncdate = yr*10000 + mon*100 + day
      call get_curr_time(ndcur, nscur)
!
!-----------------------------------------------------------------------
! Dispose history files.
!-----------------------------------------------------------------------
!
! Begin loop over ptapes (the no. of declared history files - primary
! and auxiliary).  This loop disposes a history file to Mass Store
! when appropriate.
!
      do t=1,ptapes
         if (nflds(t) == 0) cycle

         hdispose(t) = .false.
         lfill(t) = .false.
!
! Find out if file is full
!
         if (hstwr(t) .and. nfils(t) >= mfilt(t)) then
            lfill(t) = .true.
         endif
!
! Dispose history file if 
!    1) file is filled or 
!    2) this is the end of run and file has data on it or
!    3) restarts are being put out and history file has data on it
!
         if (lfill(t) .or. (nlend .and. nfils(t) >= 1) .or. (rstwr .and. nfils(t) >= 1)) then
!
! Dispose history file
!
            hdispose(t) = .true.
!
! Is this the 0 timestep data of a monthly run?
! If so, just close primary unit do not dispose. 
!
            if (masterproc) write(iulog,*)'WRAPUP: nf_close(',t,')=',trim(nhfil(t))
            if(pio_file_is_open(tape(t)%File)) then
               if (nlend .or. lfill(t)) then
                  do f=1,nflds(t)
                     if (associated(tape(t)%hlist(f)%varid)) then
                        deallocate(tape(t)%hlist(f)%varid)
                        nullify(tape(t)%hlist(f)%varid)
                     end if
                  end do
               end if
               call PIO_CloseFile (tape(t)%File)
            end if
            if (nhtfrq(t) /= 0 .or. nstep > 0) then

! 
! Print information concerning model output.
! Model day number = iteration number of history file data * delta-t / (seconds per day)
! 
               tday = ndcur + nscur/86400._r8
               if(masterproc) then
                  if (t==1) then
                     write(iulog,*)'   Primary history file'
                  else
                     write(iulog,*)'   Auxiliary history file number ', t-1
                  end if
                  write(iulog,9003)nstep,nfils(t),tday
                  write(iulog,9004)
               end if
                  !                      
! Auxilary files may have been closed and saved off without being full. 
! We must reopen the files and position them for more data.
! Must position auxiliary files if not full
!              
               if (.not.nlend .and. .not.lfill(t)) then
                  call cam_PIO_openfile (tape(t)%File, nhfil(t), PIO_WRITE)
                  call h_inquire(t)
               end if
            endif                 ! if 0 timestep of montly run****
         end if                      ! if time dispose history fiels***   
      end do                         ! do ptapes
!
! Reset number of files on each history tape
!
      do t=1,ptapes
         if (nflds(t) == 0) cycle
         lhfill = hstwr(t) .and. nfils(t) >= mfilt(t)
         lhdisp = lhfill .or. (nlend .and. nfils(t) >= 1) .or. &
              (rstwr .and. nfils(t) >= 1)
         if (lhfill.and.lhdisp) then
            nfils(t) = 0
         endif
      end do
      return
9003  format('    Output at NSTEP     = ',i10,/, &
             '    Number of time samples on this file = ',i10,/, &
             '    Model Day           = ',f10.2)
9004  format('---------------------------------------')
   end subroutine wrapup


   integer function gen_hash_key(string)
!dir$ INLINENEVER gen_hash_key
!
!-----------------------------------------------------------------------
!
! Purpose: Generate a hash key on the interval [0 .. tbl_hash_pri_sz-1]
!          given a character string.
!
! Algorithm is a variant of perl's internal hashing function.
!
!-----------------------------------------------------------------------
!
   implicit none
!
!  Arguments:
!
   character(len=*), intent(in) :: string
!
!  Local.
!
   integer :: hash
   integer :: i

   hash = gen_hash_key_offset

   if ( len(string) /= 19 ) then
!
!     Process arbitrary string length.
!
      do i = 1, len(string)
         hash = ieor(hash , (ichar(string(i:i)) * tbl_gen_hash_key(iand(i-1,tbl_max_idx))))
      end do
   else
!
!     Special case string length = 19
!
      hash = ieor(hash , ichar(string(1:1))   * 61)
      hash = ieor(hash , ichar(string(2:2))   * 59)
      hash = ieor(hash , ichar(string(3:3))   * 53)
      hash = ieor(hash , ichar(string(4:4))   * 47)
      hash = ieor(hash , ichar(string(5:5))   * 43)
      hash = ieor(hash , ichar(string(6:6))   * 41)
      hash = ieor(hash , ichar(string(7:7))   * 37)
      hash = ieor(hash , ichar(string(8:8))   * 31)
      hash = ieor(hash , ichar(string(9:9))   * 29)
      hash = ieor(hash , ichar(string(10:10)) * 23)
      hash = ieor(hash , ichar(string(11:11)) * 17)
      hash = ieor(hash , ichar(string(12:12)) * 13)
      hash = ieor(hash , ichar(string(13:13)) * 11)
      hash = ieor(hash , ichar(string(14:14)) * 7)
      hash = ieor(hash , ichar(string(15:15)) * 3)
      hash = ieor(hash , ichar(string(16:16)) * 1)
      hash = ieor(hash , ichar(string(17:17)) * 61)
      hash = ieor(hash , ichar(string(18:18)) * 59)
      hash = ieor(hash , ichar(string(19:19)) * 53)
   end if

   gen_hash_key = iand(hash, tbl_hash_pri_sz-1)

   return

   end function gen_hash_key

!#######################################################################

   integer function get_masterlist_indx(fldname)
!
!-----------------------------------------------------------------------
!
! Purpose: Return the the index of the field's name on the master file list.
!
!          If the field is not found on the masterlist, return -1.
!
!-----------------------------------------------------------------------
!
!  Arguments:
!
   character(len=*), intent(in) :: fldname
!
!  Local.
!
   integer :: hash_key
   integer :: ff
   integer :: ii
   integer :: io   ! Index of overflow chain in overflow table
   integer :: in   ! Number of entries on overflow chain

   hash_key = gen_hash_key(fldname)
   ff = tbl_hash_pri(hash_key)
   if ( ff < 0 ) then
      io = abs(ff)
      in = tbl_hash_oflow(io)
      do ii = 1, in
         ff = tbl_hash_oflow(io+ii)
         if ( masterlist(ff)%thisentry%field%name == fldname ) exit
      end do
   end if

   if (ff == 0) then
      ! fldname generated a hash key that doesn't have an entry in tbl_hash_pri.
      ! This means that fldname isn't in the masterlist
      call endrun ('GET_MASTERLIST_INDX: attemping to output field '//fldname//' not on master list')
   end if
   
   if (associated(masterlist(ff)%thisentry) .and. masterlist(ff)%thisentry%field%name /= fldname ) then
      call endrun ('GET_MASTERLIST_INDX: error finding field '//fldname//' on master list')
   end if

   get_masterlist_indx = ff
   return
   end function get_masterlist_indx
!#######################################################################

   subroutine bld_outfld_hash_tbls()
!
!-----------------------------------------------------------------------
!
! Purpose: Build primary and overflow hash tables for outfld processing.
!
! Steps:
!  1) Foreach field on masterlist, find all collisions.
!  2) Given the number of collisions, verify overflow table has sufficient
!     space.
!  3) Build primary and overflow indices.
!
!-----------------------------------------------------------------------
!
!  Local.
!
   integer :: ff
   integer :: ii
   integer :: itemp
   integer :: ncollisions
   integer :: hash_key
   type(master_entry), pointer :: listentry
!
!  1) Find all collisions.
!
   tbl_hash_pri = 0

   ff=0
   allocate(masterlist(nfmaster))
   listentry=>masterlinkedlist
   do while(associated(listentry))
      ff=ff+1
      masterlist(ff)%thisentry=>listentry
      listentry=>listentry%next_entry
   end do
   if(ff /= nfmaster) then
      write(iulog,*) 'nfmaster = ',nfmaster, ' ff=',ff
      call endrun('mismatch in expected size of nfmaster')
   end if

   
   do ff = 1, nfmaster
      hash_key = gen_hash_key(masterlist(ff)%thisentry%field%name)
      tbl_hash_pri(hash_key) = tbl_hash_pri(hash_key) + 1
   end do

!
!  2) Count number of collisions and define start of a individual
!     collision's chain in overflow table. A collision is defined to be any
!     location in tbl_hash_pri that has a value > 1.
!
   ncollisions = 0
   do ii = 0, tbl_hash_pri_sz-1
      if ( tbl_hash_pri(ii) > 1 ) then  ! Define start of chain in O.F. table
         itemp = tbl_hash_pri(ii)
         tbl_hash_pri(ii) = -(ncollisions + 1)
         ncollisions = ncollisions + itemp + 1
      end if
   end do

   if ( ncollisions > tbl_hash_oflow_sz ) then
      write(iulog,*) 'BLD_OUTFLD_HASH_TBLS: ncollisions > tbl_hash_oflow_sz', &
              ncollisions, tbl_hash_oflow_sz
      call endrun()
   end if

!
!  3) Build primary and overflow tables.
!     i - set collisions in tbl_hash_pri to point to their respective
!         chain in the overflow table.
!
   tbl_hash_oflow = 0

   do ff = 1, nfmaster
      hash_key = gen_hash_key(masterlist(ff)%thisentry%field%name)
      if ( tbl_hash_pri(hash_key) < 0 ) then
         ii = abs(tbl_hash_pri(hash_key))
         tbl_hash_oflow(ii) = tbl_hash_oflow(ii) + 1
         tbl_hash_oflow(ii+tbl_hash_oflow(ii)) = ff
      else
         tbl_hash_pri(hash_key) = ff
      end if
   end do

!
!  Dump out primary and overflow hashing tables.
!
!   if ( masterproc ) then
!      do ii = 0, tbl_hash_pri_sz-1
!         if ( tbl_hash_pri(ii) /= 0 ) write(iulog,666) 'tbl_hash_pri', ii, tbl_hash_pri(ii)
!      end do
!
!      do ii = 1, tbl_hash_oflow_sz
!         if ( tbl_hash_oflow(ii) /= 0 ) write(iulog,666) 'tbl_hash_oflow', ii, tbl_hash_oflow(ii)
!      end do
!
!      itemp = 0
!      ii = 1
!      do 
!         if ( tbl_hash_oflow(ii) == 0 ) exit
!         itemp = itemp + 1
!         write(iulog,*) 'Overflow chain ', itemp, ' has ', tbl_hash_oflow(ii), ' entries:'
!         do ff = 1, tbl_hash_oflow(ii)  ! dump out colliding names on this chain
!            write(iulog,*) '     ', ff, ' = ', tbl_hash_oflow(ii+ff), &
!                       ' ', masterlist(tbl_hash_oflow(ii+ff))%thisentry%field%name
!         end do
!         ii = ii + tbl_hash_oflow(ii) +1 !advance pointer to start of next chain
!      end do
!   end if

   return
666 format(1x, a, '(', i4, ')', 1x, i6)

   end subroutine bld_outfld_hash_tbls

!#######################################################################

   subroutine bld_htapefld_indices
!
!-----------------------------------------------------------------------
!
! Purpose: Set history tape field indicies in masterlist for each
!          field defined on every tape.
!
! Note: because of restart processing, the actflag field is cleared and
!       then set only for active output fields on the different history
!       tapes.
!
!-----------------------------------------------------------------------
!
!  Arguments:
!

!
!  Local.
!
   integer :: f
   integer :: t

!
!  Initialize htapeindx to an invalid value.
!
   type(master_entry), pointer :: listentry

   do t = 1, ptapes
      do f = 1, nflds(t)
         listentry => get_entry_by_name(masterlinkedlist, tape(t)%hlist(f)%field%name)
         if(.not.associated(listentry)) then
            write(iulog,*) 'BLD_HTAPEFLD_INDICES: something wrong, field not found on masterlist'
            write(iulog,*) 'BLD_HTAPEFLD_INDICES: t, f, ff = ', t, f
            write(iulog,*) 'BLD_HTAPEFLD_INDICES: tape%name = ', tape(t)%hlist(f)%field%name
            call endrun
         end if
         listentry%act_sometape = .true.
         listentry%actflag(t) = .true.
         listentry%htapeindx(t) = f
      end do
    end do
            
   return
   end subroutine bld_htapefld_indices

!#######################################################################

   logical function hist_fld_active (fname)
!
!------------------------------------------------------------------------
!
! Purpose: determine if a field is active on any history file
!
!------------------------------------------------------------------------
!
! Arguments
!
      character(len=*), intent(in) :: fname ! Field name
!
! Local variables
!
      character*(max_fieldname_len) :: fname_loc  ! max-char equivalent of fname
      integer :: ff                  ! masterlist index pointer
!-----------------------------------------------------------------------

      fname_loc = fname
      ff = get_masterlist_indx(fname_loc)
      if ( ff < 0 ) then
         hist_fld_active = .false.
      else 
         hist_fld_active = masterlist(ff)%thisentry%act_sometape
      end if

   end function hist_fld_active

end module cam_history
